Implementing implicit method for system of ODEs

I am trying to implement implicit Runge Kutta for solving a system of ODEs. Currently I was able to do it for scalar ODE (y'=-y), but I got stuck on how to modify it so that it work for a system of ODEs. Something like y1' = y2, y2' = -y1.

Everythin before the loop is fine if we replace the function handle f by the following function

function dydt = f(t,y)
dydt(1) = y(2);
dydt(2) = -y(1);
end

and modify the initial value to be something like

y0 = [1;0]

But How do we modify the process of solving the nonlinear system inside the loop? Also I feel I have hard-coded the method for this specific 2 stage method. Is there any way to generalize the loop for any implicit runge kutta method with size s by s?

Here is the code

n = 50; % number of time steps
t0 = 0; % initial time
T = 10; % final time
dt = (T-t0)/n; % size of timestep
y0 = 1; % initial data
m = length(y0); % dim. of problem
% intitalize time and solution vector
t = zeros(n+1,1);
y = zeros(n+1,m);
% initial conditions
t(1,1) = t0;
y(1,:) = y0(:);
% RHS of ODE
f = @(t,y) -1 * y;
% Runge Kutta coefficients
A = [1/2 + 1/(2*sqrt(3)) , 0;
     -1/sqrt(3), 1/2 + 1/(2*sqrt(3))];
b = [1/2;1/2];
c = [1/2 + 1/(2*sqrt(3));1/2 - 1/(2*sqrt(3))];
s = length(c);
for i = 1:n  

    XI = @(xi) [xi(1) - (y(i,:) + dt*A(1,1)*f(t(i)+c(1)*dt,xi(1))...
                                + dt*A(1,2)*f(t(i)+c(2)*dt,xi(2)));
                xi(2) - (y(i,:) + dt*A(2,1)*f(t(i)+c(1)*dt,xi(1))...
                                + dt*A(2,2)*f(t(i)+c(2)*dt,xi(2)))];
    xi = fsolve(XI,zeros(s,1));
    t(i+1,1) = t(i,1) + dt; 
    y(i+1,:) = y(i,:) + dt * (b(1)*f(t(i)+c(1)*dt,xi(1))+...
                              b(2)*f(t(i)+c(2)*dt,xi(2)));
end
plot(t,y)
xlabel('t')
ylabel('y(t)')

Answer 1

[Torsten]

n = 50; % number of time steps
t0 = 0; % initial time
T = 2*pi; % final time
dt = (T-t0)/n; % size of timestep
y0 = [1;0]; % initial data
m = length(y0); % dim. of problem
% intitalize time and solution vector
t = zeros(1,n+1);
y = zeros(m,n+1);
% initial conditions
t(1) = t0;
y(:,1) = y0;
% RHS of ODE
f = @(t,y) [y(2);-y(1)];
% Runge Kutta coefficients
A = [1/2 + 1/(2*sqrt(3)) , 0;
    -1/sqrt(3), 1/2 + 1/(2*sqrt(3))];
b = [1/2;1/2];
c = [1/2 + 1/(2*sqrt(3));1/2 - 1/(2*sqrt(3))];
s = length(c);
xiguess = zeros(m*s,1);
for i = 1:n  
    XI = @(xi) [xi(1:m)     - f(t(i)+c(1)*dt,y(:,i)+dt*(A(1,1)*xi(1:m)+A(1,2)*xi(m+1:m*s)));   
                xi(m+1:m*s) - f(t(i)+c(2)*dt,y(:,i)+dt*(A(2,1)*xi(1:m)+A(2,2)*xi(m+1:m*s)))];
    xi = fsolve(XI,xiguess);
    t(i+1) = t(i) + dt; 
    y(:,i+1) = y(:,i) + dt * (b(1)*xi(1:m) + b(2)*xi(m+1:m*s)); 
    xiguess = xi;
end

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

plot(t,y)
xlabel('t')
ylabel('y(t)')