Matlab - Lotka-Volterra

function [t,u,v]=forward_euler(f1, f2,u0,v0, t0, tf, n)

if nargin<1

   syms x y

    alpha=1.5; betha=2.8;

    gama=3.4; delta=1.9;

    f1=alpha*x-betha*x*y;

    f2=-gama*x+delta*x*y;

    tf=2;

    t0=0.2;

    n=10;

    u0=0.2;

    v0=0.1;

end

fun1=str2func(strcat('@(x,y,t)',char(f1)));

fun2=str2func(strcat('@(x,y,t)',char(f2)));

h=(tf-t0)/n;

u(1)=u0;

v(1)=v0;

t(1)=t0;

    for i=2:n+1

        deriv1=feval(fun1,u(i-1),v(i-1),t(i-1));

        deriv2=feval(fun2,u(i-1),v(i-1),t(i-1));

        u(i)=u(i-1)+h*deriv1;

        v(i)=v(i-1)*h*deriv2;

        t(i)=t(i-1)+h;

    end

        plot(t,u,t,v)

    grid on

    return;

end

> [t,u,v]=forward_euler

t =

    0.2000    0.3800    0.5600    0.7400    0.9200    1.1000    1.2800    1.4600    1.6400    1.8200    2.0000

u =

    0.2000    0.2439    0.3112    0.3950    0.5017    0.6371    0.8091    1.0276    1.3050    1.6573    2.1048

v =

    0.1000   -0.0116    0.0017   -0.0003    0.0001   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000

>> [t,u,v]=forward_euler('2*x+y^2+t','5*x+4*y+t/2',0.5, 0.8, 0.1, 3, 10)

t =

    0.1000    0.3900    0.6800    0.9700    1.2600    1.5500    1.8400    2.1300    2.4200    2.7100    3.0000

u =

  1.0e+215 *

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    3.7772

v =

  1.0e+216 *

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    1.8886

function [t,u,v]=forward_euler(f1, f2,u0,v0, t0, tf, n)

if nargin<1

    f1='1.5*x-2.8*x*y';

    f2='-3.4*x+2.8*x*y';

    tf=2;

    t0=0.2;

    n=10;

    u0=0.2;

    v0=0.1;

end

fun1=str2func(strcat('@(x,y,t)',char(f1)));

fun2=str2func(strcat('@(x,y,t)',char(f2)));

h=(tf-t0)/n;

u(1)=u0;

v(1)=v0;

t(1)=t0;

    for i=2:n+1

        deriv1=feval(fun1,u(i-1),v(i-1),t(i-1));

        deriv2=feval(fun2,u(i-1),v(i-1),t(i-1));

        u(i)=u(i-1)+h*deriv1;

        v(i)=v(i-1)*h*deriv2;

        t(i)=t(i-1)+h;

    end

        plot(t,u,t,v)

    grid on

    return;

end

>>  [t,u,v]=forward_euler

t =

   0.20000   0.38000   0.56000   0.74000   0.92000   1.10000   1.28000   1.46000   1.64000   1.82000   2.00000

u =

   0.20000   0.24392   0.31116   0.39491   0.50160   0.63701   0.80901   1.02743   1.30484   1.65715   2.10458

v =

  1.0000e-001  -1.1232e-002  1.6922e-003  -3.2180e-004  7.7794e-005  -2.3879e-005  9.3095e-006  -4.6092e-006  2.8982e-006  -2.3144e-006  2.3472e-006

>>  [t,u,v]=forward_euler('2*x+y^2+t','5*x+4*y+t/2',0.5, 0.8, 0.1, 3, 10)

t =

   0.10000   0.39000   0.68000   0.97000   1.26000   1.55000   1.84000   2.13000   2.42000   2.71000   3.00000

u =

   5.0000e-001   1.0046e+000   2.2164e+000   8.5335e+000   3.2694e+002   8.0343e+005   4.6418e+012   1.5646e+026   1.7743e+053   2.2826e+107   3.7772e+215