@Preben Alsholm 

Thank you for your remarks.

I have some question in the same code.

Here the code

> restart;
> with(PDEtools);
> with(LinearAlgebra);
> ode := diff(u(t), t, t)+2*GAMMA*(diff(u(t), t))+omega^2*u(t) = 0;
> deq1 := diff(u(t), t) = v(t);
> deq2 := subs(deq1, ode);
> dsolve({deq1, deq2}, {u(t), v(t)});
> eqns := [rhs(deq1) = lhs(deq1), rhs(deq2) = lhs(deq2)];
> y := [u, v]; b := diff(y(t), t);
> A, b := GenerateMatrix(eqns, y(t));
> gnat := Eigenvectors(A);
> lambda := gnat[1]; Lambda := gnat[2];
> Y := Vector([y]);
> x := [x1, x2];
> tr := GenerateEquations(Lambda, [x[1], x[2]], Y);
> tr1 := solve(%, {u, v});
> tr := subs(x[1] = x[1](t), x[2] = x[2](t), u = u(t), v = v(t), tr1);
> dchange(tr, {deq1, deq2}, [x[1], x[2]]);
> sys := solve(%, {diff(x[1](t), t), diff(x[2](t), t)});
> sol := dsolve(sys);
### In All the next I would to put Gamma=0; But How!!!!!
> x2_sol := unapply(rhs(sol[2]), t);
> x1_sol := unapply(rhs(sol[1]), t);
> f := proc (t, x, alpha) options operator, arrow; alpha*x end proc;
> Stencil[1] := x_new-x_old-f(t_old, x_old, lambda[1])*h = 0;
> Stencil[2] := x_new-x_old-f(t_old, x_old, lambda[2])*h = 0;
> stencil := convert(Stencil, list);
> Labels := [`Forward Euler `, `Forward Euler`];
> for i to 2 do

Labels[i], stencil[i]

end do;
>
> ### We want to convert these to numerical algorithms, so we need to write these as mappings:

> Sol := proc (omega, h, T, choice)

local N, t, x, i;

N := round(T/h)+1;

t[0] := 0;

x[0] := 1;

for i to N do

t[i] := t[i-1]+h;

x[i] := stencil[choice](h, omega, x[i-1])

end do;

[seq([t[i], x[i]], i = 0 .. N)]

end proc;
>
## I would like to give a condition on the eignevalue of Lambda to get the stability.
> # I would like to give a condition and plot the region where Th Euler Method is stable.

>
>
>

Thanks.

@Carl Love 

Ok. I agree. Thanks

@sarra 

Hi. Thanks for your remarks.

can we add in the previous code,  stentiel:=y(x+h)= series (y(x+h),h=0,N),

@sarra 

Thanks for your remarks, and your idea.

I try using another method,  please how to simplify this quantities in maple.

I used   simplify(eq7), but there is no simplification I used olso lhs(eq7)=simplify(rhs(eq7)), the same conclusion.

@Carl Love 

Thanks for your  answer. I have a simple, question, is it possible to know that the order of the  rest is exacty 3, using Maple.

@Joe Riel 

For me I think That F(x)= sum_{-infty}^{infty} Cn(f) exp(I*n*x)

avec Cn(f)=  (I/n-1)*(-1)^{n-1}  if n different de 1  If n=1, c1(f)=0

@Preben Alsholm 

Thanks you very much for the code and the link.

I tried  it It's working fine.

@sarra 

Dear all,

To compute the adams Moultin method, I get two Method. Which from them correspond to my question,  and I don't understant the procedure proposed. 

the First Method:

> Adamsmoulton := proc (k::posint) local P, t, f, y, n;

P := interp([t[n]+h, seq(t[n]-j*h, j = 0 .. k-2)], [f[n+1], seq(f[n-j], j = 0 .. k-2)], x);y[n+1] = y[n]+simplify(int(P, x = t[n] .. t[n]+h))

end proc;

Second Method:

@Markiyan Hirnyk 

In the link, they solves IVP using

but for me, i want a    proc(f, n,N)!!!!

I have not any idea.

Thanks

@Preben Alsholm 

Sorry, when I said it was my last question, for me it's the last question in the  another suject. I get the full answer.   Again, you have helped me and thank you very much.

@sarra 

Maye the last question, thanks Preben ALsholm.

When I put this in my code in Maple:

> taylor(eval(rhs(res), x = h)-rk3, h = 0);

I get  (1/24)*h^4-(1/120)*h^5+O(h^6)

The question is : RK3 is a method of order O(h^4) , but the taylor expansion of the error give: (1/24)*h^4-(1/120)*h^5+O(h^6) !!!  You think there is a mistake in my algorithm or what???

@Preben Alsholm 

Thank you very much.   The question is solved.

Thank your for your answer.

My goal is to Solve ODE
                             d                   
                            --- y(x) = f(x, y(x))
                             dx                  

Using 3-step Runge Kutta  Method.

> RK3 := proc (f, a, b, y0, N)

local x, y, n, h, k;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf(b-a)/N;

x[0] := a;

y[0] := 1;

for n from 0 to N-1 do

x[n+1] := a+(n+1)*h;

k[1] := f(x[n], y[n]);

k[2] := f(x[n]+(1/2)*h, y[n]+(1/2)*h*k[1]);

k[3] := f(x[n]+h, y[n]+h*(-k[1]+2*k[2]));

y[n+1] := y[n]+(1/6)*h*(k[1]+4*k[2]+k[3])

end do;

[seq([x[n], y[n]], n = 0 .. N)]

end proc;

An example is to solve : y'(x)=-y,  with intial condition y(0)=1

 ode := diff(y(x), x) = -y(x);

res := dsolve({ode, y(0) = 1});

we get as exact solution :

My Question, how to compute the error between exact and approximation solution with RK3 and show that the eroor is O(h^4) locally and globaly O(h^3).

@Preben Alsholm 

Thanks I agrre, I test your code, it's nice.

But, for me:

The n-step Adams Moulton method to solve y'(x)=F(x,y(x)) is defined by the stencil

y(x+h)=y(x)+h *sum_{j=-1}^{n-1} beta_j F( x-j*h, y(x-j*h)  )   + O(h^{n+2})

I want a procedure with single argument ''n'' that calculates and return the ''beta_i'' coefficients

Is it possible with your method. Thanks or any other idea.

@sarra 

First, special thanks to @Preben Alsholm  for all the remarks.
Now, I would like to writte a procedure that solves for y(1) in the initial value provlem:  y'(x)=f(y), y(0)=1,  using a numerical stencil based on the n^{eme} order taylor series expansion of y. The procedure should include a function f, an integrer n representing a accuracy of the Taylor series expansion, and N represent the number of steps between x=0 and x=1.

Thanks many thanks