Euler2 (a custom-defined procedure for two-dimensional systems)

DO NOT CHANGE ANYTHING IN THE FOLLOWING INPUT REGION - YOU HAVE BEEN WARNED !!!

> Euler2 := proc( MODEL, IC, DOMAIN, N )

> local a, b, dt, dx, dy, f, g, i, ode1, ode2, tt, xx, yy, LISTpts, vars;

> ode1 := op(1,MODEL);

> ode2 := op(2,MODEL);

> vars := (op(2,lhs(ode1) ), op([1,0],lhs(ode1) ), op([1,0],lhs(ode2) ) );

> f := unapply( rhs(ode1), vars );

> g := unapply( rhs(ode2), vars );

> a := op(1,rhs(DOMAIN)); b := op(2,rhs(DOMAIN));

> dt := (b-a)/N;

> tt := a;

> xx := subs( IC, vars[2](a) );

> yy := subs( IC, vars[3](a) );

> #yy := rhs(IC);

> LISTpts := [ tt, xx, yy ];

> for i from 1 to N do

> dx := evalf( f(tt,xx,yy) ) * dt;

> dy := evalf( g(tt,xx,yy) ) * dt;

> xx := xx + dx;

> yy := yy + dy;

> tt := tt + dt;

> LISTpts := LISTpts, [ tt, xx, yy ];

> od;

> RETURN( [ LISTpts ] );

> end:

DO NOT CHANGE ANYTHING IN THE PRECEDING INPUT REGION - YOU HAVE BEEN WARNED !!!

>

Let's recall what we are working with.

> MODEL;

> IC;

> VAR;

> DOMAIN ;

The number of steps to be used in Euler's method is

> Nsteps := 10;

>

The points generated by Euler's method are

> ptsE := Euler2( MODEL, IC, DOMAIN, Nsteps );

>

These data points can be used to create an approximate phase portrait as follows:

> ptsXY := [ seq( [ pt[2], pt[3] ], pt=ptsE ) ]:

> XYplot := plot( ptsXY, color = GREEN, thickness = 2 ):

> XYplot;

>

Another way to view the approximate solution is to overlay the computed solution on the direction field (and Maple's phase portrait).

> display( [ DIRplot, PHASEplot, XYplot ] );

>

Now, increase the time domain to show the solution for t=25. Recompute the Euler solution with 10 steps (dt=0.25). Repeat this problem with initial condition x(0)=2, x'(0)=0. What can you conjecture about all solutions? (Consider more initial conditions if necessary.)