Using bvp4c for coupled second order ODEs. Issue in calling bvp4c.

I have a system of two 2nd order ODEs. The set up of the ODE system is in the function bvp_function. I get the following error messages when I run the code:

Error using bvparguments
Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT):
The derivative function ODEFUN should return a column vector of length 100.

Error in bvp4c (line 122)
bvparguments(solver_name,ode,bc,solinit,options,varargin);

Error in NewPoiseuilleBaseFlow (line 38)
sol = bvp4c(@(y,X) bvp_function(y,X,b), @(Xa,Xb) boundary_conditons_fcn(Xa,Xb), guess);

Here is the code. Anyone know how to solve this problem?

% Laminar base flow for Poiseuille flow

% using Deka 2023

clear
clc

%% Setting variables

N = 100; % number of collocation points
Mach = 2.0; % Mach number
rho0 = 1.225; % density [kg/m^3]
h = 1; % height of half channel (full channel height = 2h) [m]
gamma = 1.4;
R = 287;
T0 = 288.5; % temperature (and wall temperature) (15 degC) [K]
mu0 = 1.802e-5; % dynamic viscosity [kg m^-1 s^-1]
nu0 = mu0/rho0; % kinematic viscosity [m^2/s]
kappa0 = 0.02476; % thermal conductivity [W/mK]
u0 = Mach*sqrt(gamma*R*T0); % speed [m/s]
Re = rho0*u0*2*h/mu0; % reynolds number
Pr = gamma*R*mu0/((gamma-1)*kappa0); % Prandtl number
b = Pr*(gamma-1)*Mach^2; % constant for setting up ode system

y = linspace(-1,1,N)'; % y vector
T = linspace(T0,T0,N)'; % T vector
u = linspace(1,1,N)'; % initial u vector

%% Solving the x momentum and energy equations using ode45

% initial guess for the solver
uinit = incompVelProfile(N,Mach)';

guess = bvpinit(y, uinit);

%% Solving the system of ODEs using bvp4c

sol = bvp4c(@(y,X) bvp_fn(y,X,b), @(Xa,Xb) BCs(Xa,Xb), guess);

%% Function setting up matrix system of ODEs

function dxdy= bvp_fn(y,X,b)
% X(1)=u, X(2)=u', X(3)=T, X(4)=T' 
    
    dxdy(1) = X(2);
    dxdy(3) = X(4);             
    dxdy(4) = -0.5*X(3)*(X(4))^2-b*(X(2))^2;
    dxdy(2) = -2*(X(3))^(-3/2)-0.5*(X(3))^(-1)*X(4)*X(2);
    dxdy = dxdy(:);
end

%% Function for boundary conditions

function res = BCs(Xa,Xb)
% u=0, T=1 at y=-1 and u=0, T=1 at y=1.
    res = [Xa(1)-0; 
           Xb(1)-0
           Xa(3)-1
           Xb(3)-1]; 
end

%% Function for incompressible velocity profile

function [u] = incompVelProfile(N,Mach)
% N = number of collocation points
% Mach = desired Mach number

% returns the incompressible velocity profile for channel flow

% setting variables

rho = 1.225; % density of air [kg/m-3]
nu = 1.45e-5; % kinematic viscosity [m^2/s]
T = 288; % temperature of air, 15 degC [K]
gamma = 1.4; % ratio of specific heat capacities
R = 287; % universal gas constant for air [J/(kg.K)]

% setting domain

y = linspace(-1,1,N); % y points

% working out the equivalent pressure gradient

c = sqrt(gamma*R*T); % speed of sound [m/s]
umax = Mach*c; % max velocity (centre of channel) [m\s]
dPdx = -2*umax*rho*nu; % pressure gradient required for centreline vel [Pa/m]

% obtaining velocity profile

u = zeros(N); % initialising vector

u = -0.5.*1./(rho.*nu).*dPdx.*(1-y.^2);

% have to normalise it
u = u./umax;

end

Answer 1

[Torsten]

guess = bvpinit(y, uinit);

y must be the spatial mesh,

uinit must either be a function handle that returns a 4x1 vector of initial values for the 4 solution components of X with an y value as input

or

a constant 4x1 vector that is assigned as initial value for the 4 solution components for each y value.

In your case, uinit seems to be an NxN matrix.

Read here about what uinit can be set to:

https://uk.mathworks.com/help/matlab/ref/bvpinit.html