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 := %?;

>

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 ] );

>

In a similar way, the two solution curves can be created:

> ptsTX := [ seq( [ pt[1], pt[2] ], pt=ptsE ) ]:

> plotTX := plot( ptsTX, color = GREEN, thickness = 2 ):

> plotTX;

>

> ptsTY := [ seq( [ pt[1], pt[3] ], pt=ptsE ) ]:

> plotTY := plot( ptsTY, color = GREEN, thickness = 2 ):

> plotTY;

>

And, again together with Maple's solution curve

> display( [ plotTX, P1 ] );

> display( [ plotTY, P2 ] );

>