Solving ODEs using trapezoidal method in vector/matrix notation

1 view (last 30 days)
Brenton Foreman on 22 Oct 2021
Answered: Chris on 22 Oct 2021
I have been given an example of solving multiple differential ODEs using the numerical trapezoidal method, and encountering issues when attempting to index vector and matrix values inside a for loop.
Here is the example: My code is as follows:
clear; clc;
Tend = 0.1; %End time
N = 1000; %Number of steps
T = Tend/N; %step size
t = 0:T:Tend; %time span
A = [0 -5.38e6; 1 -47.106]; %System Matrix
B = [5.38e6; 0]; %Input Matrix
I = eye(2);
C1 = inv(I - A.*(T/2));
C2 = (I + A.*(T/2));
C3 = C1.*C2; %x(t) coefficient
C4 = (B.*(T/2)); %u(t) coefficient
% Initial Conditions
x1(1) = 0;
x2(1) = 0;
%Create Cell Array
x = zeros(2,1);
xx = size(x);
y = cell(xx);
%Functions
xdot1 = @(x1,t) x1;
xdot2 = @(x2,t) x2;
xdot = @(x1,x2,t)[xdot1;xdot2];
%Assign Values to Cell Array
y{1} = xdot1;
y{2} = xdot2;
%Input
ut = @(t) 187.8.*cos(377.*t);
%Pre-Allocate and store results
xt1 = zeros(2,N+1);
%Run and store values
for n=1:N
xt1(n+1)= C3.*[y{1}; y{2}] + C4.*(ut(t(n)+ut(t(n+1))));
end
I am recieving the following error:
Error using vertcat
Nonscalar arrays of function handles are not allowed; use cell arrays instead.
Any help at finding my mistake would be greatly appreciated

Chris on 22 Oct 2021
A function handle works like this:
f = @(x,t) x.^2;
y = f(3,12);
the function f takes 3 and 12 as input variables and names them x and t, respectively. Since t is unused in the function, it's irrelevant and may not need to be in the function. It squares x and returns 9. equivalently, f = @(x) x.^2;
A couple things I see:
xdot1 = @(x1,t) x1;
You don't have a function x1 defined--x1 is a vector with an initial value. Calling xdot1(someX,someTime) would always return the entire vector x1 (initially 0). Similarly for xdot2.
xt1(n+1)= C3.*[y{1}; y{2}] + ...
y{1}, at this point, is a function
y{1} = @(x1,t) x1
It can't be concatenated into an array. If you give it numerical input, it will return a number, which is possible to use in an array.
xt1(n+1)= C3.*[y{1}(3,5); y{2})4,2)] + ...
Does this help? I'm sure there will be some more work to do once these errors are cleared up.