@Carl Love 

If N=500 (number of  iterations). I want to save only 10 of P[k].

In C I do :  If (k%(N/10)==0) .......... 

@Carl Love 

Thanks.

How can I do it if have to display plots of <X|A> for each time step .

@Rouben Rostamian  

Hi, Sorry to get you back to my problem.

See below my code about solving  my PDE at the first time step.

I want to plot A[j] as a function of X[j]. How to do this? Thanks in advance.

@tomleslie 

As I've shown above and for more explanation,

After computing A[] we  assign a[]=A[] and continue iteration.

At the end, it's just a way to make difference between time step  and next time step  so as not loosing the information localized in a[] during its computation at next time step.

@Rouben Rostamian  

Thanks for your consideration to my problem.

Exactly, I used A[] rather than a[]  to avoid confusion between time step and next time step.  

1D Boundary condition is A[0]=A[1] and A[M]=A[M-1] if I consider the grid points step is h=1/M and the time step is tau starting from 0 to N.

Here, I'm using classical theta scheme to compare results between Implicit, explicit and Crank-Nicolson.

@tomleslie 

Sorry to bother you. Thanks for your precious time.

1- the context is a resolution of generalized riemann problem like U_t+f(U)_x=0 where U=(u,v).

2- The solution constructed is u(x,t)=u0(z)+t*u1(z)+...  where z=x/t. u0,v0 are known and u1 should be solved according to the following ODE system  dsys:={(1-v0')*z*u1-(z*u0'+u0)*v1=z*u0^m, u1'+z*v1'+v1=0} where m<0.

ICS: U(x,0):=piecewise(x>=0,1,0)

3- The analytic resolution is what was given in the worksheet: resolution of v1 from eq1:=-z.v1'+(w1*v1)'=w2' and then extract u1.

4- What is expected : u positive, u -> infnity when z ->0 when t increases ,    plot numeric solution.

I hope that it is well explained.

Thanks in advance.

@tomleslie 

I think I should review the equation u0(z)+t*u1(z) so as to extend it in power 2 or 3 of t

until I could (if it's possible) reach positive solution.

Thanks anyway.

@tomleslie 

There was a mistake on v0. I changed to:

v0:= proc(z) piecewise(z<0,-(2/3)*z+c1/sqrt(-z),-(2/3)*z+c2/sqrt(z)); end proc:

Then, I want to subdivide my solution on three regions:

z<0 -> eq30, z>0 -> eq31 and (eqc20, eqc21) when solutions of eq30 and eq31 become negative or complex.

@Preben Alsholm 

Thanks.

But, there is something wrong in my code. plot shows the solution becoming negative in an interval.

How to extract values of z where solutions of eq1 switch to négative side so as I can integrate eqc  ?

(It's clear that z=6^(2/3)/4 and z=-6^(2/3)/4 are worong) .

I mean without picking up manually the values at each changing of parameters.

Thanks in advance.

@tomleslie 

Thanks a lot tomleslie. You've made life easier.

Thanks also to the all intervenors.

Best Regards.

@Preben Alsholm

Hi,

If you could help us to find what it is exactly rfk45_dae. Tomeslie said that it will be a  wasting time to seek a rfk45 algorithm to my problem solved using this method.

I have to write at least the principles steps that given by rfk45_dae method. It will be so difficult to provide a code without knowing what is exactly. 

Thanks for your habitual help.

@tomleslie 

Hi,

According to the code below, rfk45 method has solved the problem.

My question is how can I rewrite the principle steps that Maple has done in this case.

I thought to the conversion given in you reply. If not, what was the algorithm ?

Thanks in advance.

worksheet.mw

@tomleslie 

If you could help me to determine rfk45 steps. As It is known that at each step of RFk order 4, the final answer is computed using the formula

So, I would like to write the formula of each ki. To do this, I think that U' should not be function of a derviative parameter  U'=G(x,U,S,F,F').

I hope that I success to explain my problem.

@tomleslie 

See below my system solved using rfk4_dae method.

Thanks.

worksheet.mw

@Preben Alsholm

Sorry, see attached the worksheet.

worksheet.mw