@mehdi jafari 

Thanks for all the remarks about this question. 

## The code is very nice, now it's working. For me, I note u[i,j]^k=u(x(i),y(j),t(k));

N:=2;

x := i -> 1/(N)*i;
y := j -> 1/(N)*j;
t := k -> 1/(N)*k;
x[0] = x(0),x[N] = x(N),y[0] = y(0),y[N] = y(N),t[0] = t(0),t[N] = t(N);
#Here is a visualization of the lattice:
Zx := [seq]( x(i),i=0..N);
Zy := [seq]( y(i),i=0..N);
Zt := [seq]( t(i),i=0..N);

Tab:=[seq(seq(seq([Zx[i],Zy[j],Zt[k]],i=1..N+1),j=1..N+1),k=1..N+1)];
pointplot3d(Tab,color = magenta,axes = boxed, symbol = box);


Question: How can I put at each point  (x[i],y[j],t[k]) its name  u[i,j]^k, and  tickmarks in the x and y axis.

I try these lines but I think there is a problem:

textplot([seq([i,-0.1,typeset(u(x+i*h,y,t))],i=-1..1),[0,0.9,typeset(u(x,y+h))],[0,-1.1,typeset(u(x,y-h))]],align={below})],view=[-1.4..1.4,-1.4..1.4],tickmarks=[[-1=x-h,0=x,1=x+h],[-1=y-h,0=y,1=y+h]]

Many thinks.

@mehdi jafari 

Now, It's okay that is I need. Many thanks.

@sarra 

The seq of all point is:
x := [seq]( (1/5)*i,i=0..12);  #  x[i] the x-coordinate
     

y := [seq]( (1/5)*j,j=0..12); # y[j] the y-coordinate
      

t := [seq]( (1/5)*k,k=0..12);  #  t[k] the t-coordinate

# This table contains all the point must be ploted.

Tab:=seq(seq(seq([x[i],y[j],t[k]],i=1..3),j=1..3),k=1..3);

Now, how can we put these point in the space (x,y,t).

@mehdi jafari 

Thank you for your remark. But, it's seem for me that your plot is only the point P[i] with coordinate (x[i],y[i],t[i]).  But for me, I want all the point P[i,j,k] with coordinate (x[i],y[j],t[k]).

For example, when you put k=1, you will get all the point  P[i,j,1]  with coordinate (x[i],y[j],t[1]) for all i,j and those point is in the horizontal plane fixed by t=t[1].

I think we will get many horizontal surface caracterized by, t=t[k], and in these plane we put all the point P(i,j,k=fixed).

@Carl Love 

Thank you very much, There is a problem of stability. I just change the value of c=1/2, I have not any problem. I get the same figure.

the condition of stability is : c1^2+c2^2<=1. In my old file, I work with c1=1 and c2=1, so that the method is unstable. For example, one van take  c=1/2 , we have not any problem. 

This night  I will try to finish the code of stability.  

one again, thank you very much.

Please, if I want to plot the error ( between the exact and numerical ), how can I make the difference between the two matrix and plot the error at each time

@Carl Love 

I hope it's working now. Thank your for your help.

Analytic_numerical_s.mw

@Carl Love 

Thank you for your remark. I agree with your remark.

But, my solution is u(x,y,t), when i fixe t from Vect_T:=seq(k*0.1,k=0..10):   and i plot the two solution it's seem for some t, the two curves are the same, but if one take for exemple t=Vect_T[5], i.e working with t=0.4 the two curves are not the same also, take t=t=Vect_T[10], the same conclusion.  Please, run the following code.

Analytic_numerical_s.mw

@Preben Alsholm 

Thank you very much for the answer.  It is goal.

Many thinks.

Now, I will writte a code to compute the exact solution, then compute the error.  

@mehdi jafari 

First special thanks for your remarks.

Le me recall that our function is a function with three varibles (x,y,t)  so the curve is a surface in R^3 and not in R^2. Can we plot the curve in R^3 and Not in R^2.

u(i,j,k) is an approximation of u(x,y,t) at point (xi,yj,tk) can we plot all the point in R^3.

For example, we can fix two times, (i.e fix k)  and plot the solution for two different times.

Please see the last page in the following link:

http://pauli.uni-muenster.de/tp/fileadmin/lehre/NumMethoden/WS0910/ScriptPDE/Wave.pdf

There is two differents curves for two time in the last page in the link.

How can I have the same figure.

Thank you very much for your help.

Many thanks.

@mehdi jafari 

Thank you for your answer, I agrre with you,  I forget to put the value of c.  Recall that "u(x,y,t)" is the exact solution of my problem. After discretisation (finite difference):

u[i,j,k] is an approximation  of the analytic  solution  at point (x(i), y(j),t(k)) with

i=1..Fx; j=1..Fy, k=1..Ft.   I want to plot this Numerical  solution. 

If I use Matlab, for me I can plot the Numerical solution, for exmaple I can do:

for i= 2:Ft + 1
figure(1);
mesh(y, x, u(:,:,i));
title(‘Numerical solution’);
end

Thanks.

@Carl Love 

I make the change, but now my code run without stop, afther five minutes I Interrupte the computation, maye be a collocation memory needed, or a big error let the code run without stopping.

restart:
ode := diff(y(x), x) = 2*x-y(x);
f:=(x,y)->2*x-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));
                             d                   
                            --- y(x) = 2 x - y(x)
                             dx                  
(x, y) -> 2 x - y
                            -2 + 2 x + 3 exp(-x)
#starting and ending values of a<x<b and local error tolerance
a:=0; b:=1; epsilon:=1e-2;
## Procedure
test := proc (f, a, b, epsilon, N)
local n, x, y,k,z,R,h,hmin, hmax,hstar,p;
## Initial condition
 x[0] := 0; y[0] := 1; z[0]:=1;
  R[0]:=0; h[0]:=evalf(b-a)/N; hmax:=h[0];
hmin:=(1e-4)*hmax;


##calculate the RK steps
p:=2;
for n from 0 to N-1 do
x[n+1]:= x[n]+(n+1)*h[n];
 k[1] := f(x[n], y[n]);
 k[2] := f(x[n]+h[n], y[n]+h[n]*k[1]);
 k[3] := f(x[n]+h[n]/2, y[n]+h[n]/4*(k[1]+k[2]));
## 2-stage runge Kutta.
z[n+1] := z[n]+(h[n]/2)*(k[1]+k[2]);
y[n+1] := y[n]+(h[n]/6)*(k[1]+k[2]+4*k[3]);

while x[n]<b and h[n]>hmin do
n:=n;    #set last point to the end point
if x[n]+h[n]>b then
h[n]:=b-x[n];
end
end do;

## local Error
R[n+1]:=abs(y[n+1]-z[n+1]);
hstar[n+1]:=(epsilon/R[n+1])^(1/p);
if R[n+1]<=epsilon
then
 h[n+1]:=h[n]; # we take the same stepsize
 x[n] :=x[n]+h[n];
 z[n] :=z[n+1];
 y[n] :=y[n+1];
end if;
if epsilon <=R[n+1] then h[n+1]= min(0.9*hstar[n+1],hmax); end if;
if x[n]<b  then
printf("Integration did not reach end point requested");
  end if;
end do;
#return [seq([x[i], z[i]], i = 0 .. N)];
#return[seq([x[i], y[i]], i = 0 .. N)];
return[seq([x[i], R[i]], i = 0 .. N)];
end proc:
                                 
#test(f,0,1,1e-2,4);
test(f,0,1,1e-3,4);
Warning,  computation interrupted

@Carl Love 

I used the same Method here, in Runge Kutta with setpsize varibale, but it's run with stop.  I think maybe must I allocate memory or something else.

This is the code.

restart:
ode := diff(y(x), x) = 2*x-y(x);
f:=(x,y)->2*x-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));
                             d                   
                            --- y(x) = 2 x - y(x)
                             dx                  
(x, y) -> 2 x - y
                            -2 + 2 x + 3 exp(-x)
#starting and ending values of a<x<b and local error tolerance
a:=0; b:=1; epsilon:=1e-2;
## Procedure
test := proc (f, a, b, epsilon, N)
local n, x, y,k,z,R,h,hmin, hmax,hstar,p;
## Initial condition
 x[0] := 0; y[0] := 1; z[0]:=1;
  R[0]:=0; h[0]:=evalf(b-a)/N; hmax:=h[0];
hmin:=(1e-4)*hmax;


##calculate the RK steps
p:=2;
for n from 0 to N-1 do
x[n+1]:= x[n]+(n+1)*h[n];
 k[1] := f(x[n], y[n]);
 k[2] := f(x[n]+h[n], y[n]+h[n]*k[1]);
 k[3] := f(x[n]+h[n]/2, y[n]+h[n]/4*(k[1]+k[2]));
## 2-stage runge Kutta.
z[n+1] := z[n]+(h[n]/2)*(k[1]+k[2]);
y[n+1] := y[n]+(h[n]/6)*(k[1]+k[2]+4*k[3]);

while x[n]<b and h[n]>hmin do
n:=n;    #set last point to the end point
if x[n]+h[n]>b then
h[n]:=b-x[n];
end
end do;

## local Error
R[n+1]:=abs(y[n+1]-z[n+1]);
hstar[n+1]:=(epsilon/R[n+1])^(1/p);
if R[n+1]<=epsilon
then
 h[n+1]:=h[n]; # we take the same stepsize
 x[n] :=x[n]+h[n];
 z[n] :=z[n+1];
 y[n] :=y[n+1];
end if;
if epsilon <=R[n+1] then h[n+1]= min(0.9*hstar[n+1],hmax); end if;
if x[n]<b  then
printf("Integration did not reach end point requested");
  end if;
end do;
#return [seq([x[i], z[i]], i = 0 .. N)];
#return[seq([x[i], y[i]], i = 0 .. N)];
return[seq([x[i], R[i]], i = 0 .. N)];
end proc:
                                      0
                                      1
                                    0.01
#test(f,0,1,1e-2,4);
test(f,0,1,1e-3,4);

@sarra 

Hi. Please  I need to correct my code.

restart:
ode := diff(y(x), x) = 2*x-y(x);
f:=(x,y)->-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));
                             d                   
                            --- y(x) = 2 x - y(x)
                             dx                  
F:(x, y) -> 2*x-y

#starting and ending values of a<x<b and local error tolerance
a:=0; b:=1; epsilon:=1e-2;
## Procedure
RKadaptivestepsize := proc (f, a, b, epsilon, N)
local n, x, y,k,z,R,h,hmin, hmax,hstar,p;
## Initial condition
 x[0] := a; y[0] := 1; z[0]:=1;
  R[0]:=0; h:=evalf(b-a)/N; hmax:=h;
hmin:=(1e-4)*hmax;

while (x[n])<b) &(h>hmin)   
 n:=n+1;    %set last point to the end point
    if x[n]+h[n]>b
        h[n]:=b-x[n];
    end
end
##calculate the RK steps
p:=2;
for n from 0 to N-1 do
x[n+1] := x[n]+(n+1)*h[n];
k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h[n], y[n]+h[n]*k[1]);
k[3] := f(x[n]+h[n]/2, y[n]+h[n]/4*(k[1]+k[2]));
## 2-stage runge Kutta.
z[n+1] := z[n]+(h[n]/2)*(k[1]+k[2]);
y[n+1] := y[n]+(h[n]/6)*(k[1]+k[2]+4*k[3]);
## local Error
R[n+1]:=abs(y[n+1]-z[n+1]);
hstar[n+1]:=sqrt(epsilon/R[n+1]);
if R[n+1]<=epsilon
then
 h[n+1]:=h[n]; # we take the same stepsize
 x[n]:=x[n+1];
 z[n]:=z[n+1];
  y[n]:=y[n+1];
end if
if R[n+1]:=0 then h[n+1]:= min(0.9*hstar[n+1],hmax);  end if
if (x[n]<b)
display('Integration did not reach end point requested');
  end if;
end do;
#return [seq([x[i], z[i]], i = 0 .. N)];
#return[seq([x[i], y[i]], i = 0 .. N)];
return[seq([x[i], R[i]], i = 0 .. N)];
end proc:
                                   
   
Error, `)` unexpected
RKadaptivestepsize(f,0,1,1e-2,10);
                   

@Carl Love 

Thank you for your remark. I agree with you.  I make the change.

aaa:=proc(x,y,f,N)
local bb,i;
for i  from 1 to N do

bb[i]:=f(x[i],y);

if bb[i] <=2  then  y:=3*x[i]  else  y:=x[i]^2

end if ;

end do:

end proc:

@Carl Love 

I work in the interval [0,1].

In my code, when I run it, The points x[n+1]=x[n]+n*h, dont gives all the nodes in my interval, the bigest nodes is O.34. Why. But If I remove the condition if ( conditon about R<=epsilon) I get all the nodes. Please see this problem.