Solve ODE via Midpoint rule nonlinear system

Hey ,

I have the following nonlinear system.

And I'd like to approximate the solution via the midpoint rule , that is

I choose the stepsize h=0.01. x_{n+1} denotes the solution at time t=h(n+1).

I did the following:

fun=@root2d;
x0=[20,20];
sol=fsolve(@(x)root2d(x,20,20),x0)
      
      function F=root2d(x,xold,yold)

F(1)=x(1)-xold-dt*((1/2)*(x(1)+xold)-A*(1/2)*(x(1)+xold)*(x(2)+yold));
F(2)=x(2)-yold-dt*(-(1/2)*(x(2)+yold)+B*(1/2)*(x(1)+xold)*(x(2)+yold));
end 

but I keep getting the error code

Caused by:
    Failure in initial objective function evaluation. FSOLVE cannot continue.

Actually , I was quite confused anyway when I tried to implement it since I was not sure how to deal with the initial values and it seems like i incoorporated them twice..I hope somebody can help me !

Thanks.

Answer 1

[Alan Stevens]

You can do it this way:

A = 0.01;
B = 0.01;
dt = 0.01;
t = 0:dt:20;

x = zeros(size(t));
y = zeros(size(t));

x(1) = 20;
y(1) = 20;

for i = 1:length(t)-1
    
    xnew = x(i); ynew = y(i);
    err = 1;
    while err > 10^-6
          xold = xnew; yold = ynew;
          xav = (xnew+x(i))/2;
          yav = (ynew+y(i))/2;
          xnew = x(i) + dt*(xav - A*xav*yav);
          ynew = y(i) + dt*(-yav + B*xav*yav);
          err = max(abs(xnew-xold), abs(ynew-yold));      
    end
    x(i+1) = xnew;
    y(i+1) = ynew;
    
end

plot(t,x,t,y)
xlabel('t'),ylabel('x and y')
legend('x','y')

to get:

[Mathematica]

Thank you so much !! Appreciate it a lot!