@sand15 

thanks for reply.

Remark 1: you are right. the x derivative is also considered in this transformation as diff(..,x)= -ik. so most of the k apreas in the equation comes from derivatives wrt x. So the equation in k space that I shared with you is the Fourier transform of the main equation in real space.

remark 2: in the boundary condition I wrote, the initial shape in real space is presented and this also should transform to k space.

I believe the strategy in numerical approach should be:

1 - generating a list/array of k values and corresponding x values.

2- transform initial condition to k space

3- solve equation in maple numerically in k space. so we have eta(k,t)

4- inverse transform of eta(k,t) to get zeta(x,t)

those steps I think works. but I could not apply them in Maple

I would appreciate it if you can help

@mmcdara 

I attached the full equations and relations step by step in a Maple file. i also tried to explain them. i hope it would be helpful to understand the relation and the way to solve it

file
coefficients.mw

@mmcdara

first thanks for your reply. sample data is like this

k(1/2)                          omega(1/2)

q_{12}^{(ab)2} &=& \kappa_1^{(a)2}+\kappa_2^{(b)2}+2k_1k_2\;.

  &&\zeta^{[2]}=\sum_a\sum_b\int_{-\infty}^\infty dk_1
  \int_{-\infty}^\infty dk_2\;C^{(a)}_1C^{(b)}_2\;\exp\big\{i (k_1+k_2) x-i\left[\Omega_1^{(a)}+\Omega_2^{(b)}\right] t\big\}\nonumber\\
  &&\times\bigg\{s^{(ab)}-\frac{k_1+k_2}{\Omega_1^{(a)}+\Omega_2^{(b)}}\bigg(\cosh\left[(k_1+k_2)h\right]\epsilon_{c}^{(ab)}
  +\sinh\left[(k_1+k_2)h\right]\epsilon_{s}^{(ab)} \nonumber\\
  &&+\cosh\left[q_{12}^{(ab)}h\right]\mu_c^{(ab)}+\sinh\left[q_{12}^{(ab)}h\right]\mu_s^{(ab)}\bigg)\bigg\}\label{sz2e2}
\end{align}

Actually all coefficients would be constant and one can consider a constant value. but for me the method of solving and plot them as function of x and t is important

and i also dont know how to write and implement them in Maple.

I would appreciate it if you could help me

@Rouben Rostamian  Thank you so much. It is very helpful.

@Scot Gould first of all thank you for your comment. Actually, the whole story is that there is a non-linear PDE in real space. I transformeed it in Fourier space. as a result, I have u(k,t) in Fourier space and the PDE is

$$\hat{\u_{tt}}+i*\gamma\hat{\u_t}-\omega_0^2\hat{\u_}+\frac{3}{4}k \omega_0 \hat{\u_}*\hat{\u_}= 0$$

where \hat{\u_}*\hat{\u_} means convolution. assume that the shape of the u(x,t) has at t=0 is sech(x)^2. I want to solve this equation to in fourier space and then transform the result to real space to see the evolution of u(x,t) over time. above, I just wrote the linear part of the equation. ( Actually I have no idea how to wrote the last term :). i would appreciate if you can help)
in this case the term (i*\hat{\u_t}) in Fourier space is corresponding to u_{x,t} (I mean diff(u(x,t),x,t)
       

you can see that the initial hump decays over time

I hope I could clarify the question above

@tomleslie 

Download MAPLE_sheet.mw

@tomleslie if you mean

sol := dsolve( map(op, [ode_system, ic_system]), [{u[1](t), u[2](t), u[3](t), u[4](t), u[5](t)}, numeric, output = listprocedure] );

Error, (in dsolve) too many arguments; some or all of the following are wrong: [{u[1](t), u[2](t), u[3](t), u[4](t), u[5](t)}, [{u[1](t), u[2](t), u[3](t), u[4](t), u[5](t)}, numeric, output = listprocedure]]

Dear @dharr yes, you are right.

I corrected this mistake and got another error

Error, (in dsolve) too many arguments; some or all of the following are wrong: [{u[1](t), u[2](t), u[3](t), u[4](t), u[5](t)}, numerical, output = listprocedure]

Dear @dharr, first of all, I apology for confusion.

I want to solve this differential equation

dQ/dK = -(diff(Deq, K))/(diff(Deq, Q)

Where Q is complex and K is real by default. Because of singularities, I consider K as a complex variable. therefore I can move on contour on upper half plane of complex plane. (based on your comment, it might be a contour integration). As a result, we have K and Q. then calculate omega function as

omega = -I*sqrt(pp*9.81*1e-3^3)/1e-3^2*(K^2-Q(K)^2)

and plot is the imaginary part of omega in terms of K.

thanks for taking time on the issue

@dharr Actually I want the same plot from the Deq. this is my goal. that plot is calculated in C. In maple I could not find appropriate method and configurations to get the same plot. what I said about the singularities was just for clarification about the equation and the plot. In addition in many methods I used, I get errors or warnings about singularitieslike this

"cannot evaluate the solution further right of .62444381, probably a singularity"

thank you so much

Dear @dharr 
Stiff method works. But it make a semicircle around the first singularity which is located before K=1, and don't find or consider the second one which is around 1.8
How can I modify the stiff and the code

@dharr Thanks a lot for taking time

@dharr Thanks again for taking time.

Actually the main ODE is

dDeq/dQ * dQ/dK+ dDeq/dK =0

so

dQ/dK = -(diff(Deq, K))/(diff(Deq, Q)

in my image, there are two branches. in both cases, i believe that 0/0 occures. for the first one around K=0.5 the reason is that K and Q are small and it is possible to find it by taylor expansion. in the second one where the K is around 1.9 the numenator and denominator goes to zero. in this case I think that rediual theory work. so I need complex K instead of real one. something like this images (I found them in Google)

As shown in a figure below, the taylor expansion helps when K is small. but Maple could not handle the second branch and the plot sharply increased.

I plot this in maple.

@dharr Thaks for sharing your idea and the file. it is very helpful. However, I think sth is wrong. Actually, I expected to have plots like this

where red and blue lines related to

dsolve({ode,Q(0)=1e-15+1e-15*I}, numeric, method=rkf45, output = listprocedure)
dsolve({ode,Q(0)=1e-5+Pi/2*I}, numeric, method=rkf45, output = listprocedure)

Note that the code above does not work. the image is what I have from C code. however what I have to have is that plot

but in your version, the behaviur is different

sol1 := dsolve({eval(ode2,a=1.e-5*I),Q2(0)=1.e-8+1.e-8*I}, numeric, method=rkf45, output = listprocedure);

plt_im1:= pointplot({seq([rhs(sol1(K)[1]),-Im(-I*sqrt(pp*9.81*1e-3^3)/1e-3^2*(rhs(sol1(K)[1])^2-rhs(sol1(K)[2])^2))],K=0..10,0.1)},style=line,color=blue):
display([plt_im1],size=[600,600],view=[0..10,0..500]);

I would appreciate it if you could modify the code and learn me more.

@dharr In this case there are singularities. to avoid this, I need a path in complex plan. so K is no longer real. there is also cases where 0/0 occure. in these cases K should be complex