Unable to solve differential equation with finite difference method

I am trying to solve a differential equation of the following form:

J = m * L^2 / 3. m and L are made-up values for someones weight and length of their leg. F_m (t) = 0. Because \phi is a small angle, sin(\phi(t)) can be rewritten as \phi(t)

I have the following code:

b = 0.02;
m = 80;
L = 0.7;
g = 9.81;
J = m * L^2/3;
a = 0; bet = 2;
alpha = 0; beta = 3/2;
N = 200;
p_function = @(t) -b / J * ones(1, N); 
q_function = @(t) zeros(1, N);
r_function = @(t) (-(3 * 3 * g) / (2 * m * L^3));

% findif_ODE is the function to solve using finite difference method
[tSol, YSol] = findif_ODE([a, bet], [alpha, beta], p_function, q_function, r_function, N);

%{ 
Approximation of the solution of the linear limit problem
y''(x) = p(x) y'(x) + q(x) - r(x)     a < x < b
y(a) = alpha     y(b) = beta
with the linear finite difference scheme.
%
INPUT:
I = [a, b]: extremes of the integration interval
Ybc = [alpha, beta]: boundary conditions
p(x), q(x), r(x): known terms of the differential problem (function)
N: number of internal nodes
%
OUTPUT:
Xi(1:N+2): node vector (column vector)
Yi(1:N+2): vector of approximations (column vector)
%}

function [Xi, Yi] = findif_ODE(I, Ybc, p, q, r, N)
    % Discretization step
    h = (I(2)-I(1)) / (N+1);

    % Node vector (row vector)
    Xi = linspace(I(1),I(2),N+2);
    
    %Tried to visualize the sizes to figure out the problem
    disp(['Size of Xi ' num2str(size(Xi))])
    disp(['Size of p(Xi(3:N+1)): ' num2str(size(p(Xi)))])
    disp(['Size of q(Xi(2:N+1)): ' num2str(size(q(Xi)))])
    disp(['Size of r(Xi(2:N+1)): ' num2str(size(r(Xi)))])

    % Corfficient matrix of the linear system
    A_diag = [2 + h^2 * q(Xi(2:N+1))];         % main diagonal
    A_coinf = [- 1 - h * 0.5 * p(Xi(3:N+1))];  % inferior codiagonal
    A_cosup = [- 1 + h * 0.5 * p(Xi(2:N))];    % superior codiagonal

    %A = diag(A_diag) + diag(A_coinf,-1) + diag(A_cosup,1);
    A = spdiags([A_diag', A_coinf', A_cosup'], 0:2, N, N);

    % Vector of known terms (column vector)
    B = h^2 * r(Xi(2:N+1));
    disp(['Size of B ' num2str(size(B))])
    disp(['Size of B1 ' num2str(B(1))])
    disp(['Size of B1 ' num2str((1 + h * 0.5 * p(Xi(2))) * Ybc(1))])
    B_1_hold = B(1);
    B(1) = (1 + h * 0.5 * p(Xi(2))) * Ybc(1);
    B(1) = B_1_hold;
    B(N) = B(N) + (1 - h * 0.5 * p(Xi(N+1))) * Ybc(2);
    B = B';

    % Solution of the linear system
    Yi = A \ B;

    % Output
    Xi = Xi';
    Yi = [Ybc(1); Yi; Ybc(2)];
end

The error i get is the following:

Unable to perform assignment because the left and right

sides have a different number of elements.

Error in exercise_2_1_8>findif_ODE (line 65)

B(1) = (1 + h * 0.5 * p(Xi(2))) * Ybc(1);

Error in exercise_2_1_8 (line 19)

[tSol, YSol] = findif_ODE([a, bet], [alpha, beta], p_function, q_function, r_function, N);

I have tried to diagnose it myself, but I cannot seem to figure it out.

[Alan Stevens]

Your final boundary condition of theta = 3/2 (~ 86 degrees) is not really consistent with a small angle approximation. Since you are doing a numerical calculation there is little lost by just using sin(theta).

(This doesn't help solve your error problem I'm afraid!)

[Torsten]

Your functions p and q always return vectors of length N and r returns a scalar. This is wrong in all the below calls to the functions:

A_coinf = [- 1 - h * 0.5 * p(Xi(3:N+1))];  % inferior codiagonal
A_cosup = [- 1 + h * 0.5 * p(Xi(2:N))];    % superior codiagonal
B = h^2 * r(Xi(2:N+1));
B(1) = (1 + h * 0.5 * p(Xi(2))) * Ybc(1);
B(N) = B(N) + (1 - h * 0.5 * p(Xi(N+1))) * Ybc(2);

Answer 1

[Alan Stevens]

Should be more like this I think:

b = 0.02;
m = 80;
L = 0.7;
g = 9.81;
J = m * L^2/3;
a = 0; bet = 2;
alpha = 0; beta = 3/2;
N = 200;
p = b / J; 
r = 3*g/(2*L*J);

% Approximation of the solution of the linear limit problem
% y''(x) = -p y'(x) - r*y(x)     a < x < bet
% y(a) = alpha     y(bet) = beta
% with the linear finite difference scheme.
% %
% INPUT:
I = [a, bet]; % extremes of the integration interval
% Ybc = [alpha, beta]: boundary conditions
% p, r: known terms of the differential problem (function)
% N: number of internal nodes
% %
% OUTPUT:
% Xi(1:N+2): node vector (column vector)
% Yi(1:N+2): vector of approximations (column vector)
% %}

    % Discretization step
    h = (I(2)-I(1)) / (N+1);

    % Node vector (row vector)
    Xi = linspace(I(1),I(2),N+2)';
    

    % Coefficient matrix of the linear system
    A_diag = diag(-(2-r*h^2)*ones(1,N+2));         % main diagonal
    A_diag(1,1) = 1;  A_diag(N+2,N+2) = 1;
    A_coinf = diag((1-p*h/2)*ones(1,N+1),-1);  % inferior codiagonal
    A_cosup = diag((1+p*h/2)*ones(1,N+1),1);   % superior codiagonal

    A = A_diag + A_coinf + A_cosup;
    A(1,2) = 0;  A(N+2,N+1) = 0;
    
    % Vector of known terms (column vector)
    B = zeros(N+2,1);
    B(1) = alpha;  B(N+2) = beta;
    
    % Solution of the linear system
    Yi = A \ B;
        
    % Output
plot(Xi,Yi),grid