Why do you use ode45 for a boundary value problem ? Use bvp4c instead.
The below code doesn't give errors or warnings - at least as long as it ran under MATLAB online here.
%2) loop to calculate compute boundary layers:
%% Blasius Boundary Layer Solution and CBL Solution for Flat Plate
%% Patrick D. Sinko Uppsala University, February 2, 2023
%purpose To solve the system of ODE for the hydrodynamic and concentration
%equations for flow over a flat plate
%Bisection Method Solver code modified from: Mohamad Aslani (2023).
% Compressible Similarity Solution for Flat Plate BL
% (https://www.mathworks.com/matlabcentral/fileexchange/73587-compressible
% -similarity-solution-for-flat-plate-bl), MATLAB Central File Exchange.
% Retrieved February 5, 2023.
%%
clear all
clc
close all
%% Global initialization
nu = 0.000000696;% [m^2/s] Water kinematic Viscosity at 37C
U_inf = 0.5; %m/s
Sc = 1; % starting value
% Initializing loop variables
step=1;
end_step=2;
% vector defining the span and interval of Schimdt number to be investigated.
SC = (Sc:step:end_step);
ff = zeros(length(SC),1);
gg = zeros(length(SC),1);
eta_f = zeros(length(SC),1);
eta_g = zeros(length(SC),1);
count_j =0; %storage variable counter
while(Sc <= end_step)
%% Initial guesses for f''(0): need not to be changed for most cases...
format long;
af = -0.0;
bf = 2.0;
ag = -20.0;
bg = 20.0;
BC_f = 1; % goal for the f' boundary condition i.e. (u/Uinf) = 1
BC_g = 1; % goal for the g boundary condition i.e. (c/cinf) = 1
%% Solver specific parameters
eta_max_ode = 10; % length of eta distance
% Maximum iterations
max_iter = 1e2;
% Iteration counter
count = 0; %convergence iteration while counter
% Tolerance for boundary condition
epsilon = 1e-6;
% Value greater than epsilon to start it off
error = 1e15;
ERR = zeros(max_iter,1); %initializing here clears the previous loop's data
ERR_f = zeros(max_iter,1);%initializing here clears the previous loop's data
ERR_g = zeros(max_iter,1);%initializing here clears the previous loop's data
%% Solve the ODE using ode45
tic
while(abs(error) > epsilon && count < max_iter)
sigmaf = (af+bf)/2; % initial guess
sigmag = (ag+bg)/2;
% sigmaf = 0.332;
f0 = [0, 0, sigmaf,0,sigmag];
options = odeset('RelTol',1e-9,'AbsTol',1e-9);
% the Hydrodynamic and concentration set of 1st order ODE
[eta,f] = ode45(@(t,y)Blas_Sc(t,y,Sc),[0, eta_max_ode],f0,options);
% Perform bisection method
errorf = f(end,2) - BC_f; % error in the velocity profile
errorg = f(end,4) - BC_g; % error in the concentration profile
% error = abs(errorf) + abs(errorg) + abs(errorh)
error = abs(errorf) + abs(errorg);
% error = abs(errorg);
if(errorf > 0)
bf = sigmaf;
else
af = sigmaf;
end
if(errorg > 0)
bg = sigmag;
else
ag = sigmag;
end
ERR(count+1) = error; %tracking the convergence of the total error
ERR_f(count+1) = errorf; %tracking the convergence of the hydrodynamic error
ERR_g(count+1) = errorg; %tracking the convergence of the concentration error
% Increment the counter
count = count + 1;
end
toc
count=0; %reset the counter for next loop (next input of Sc)
count_j = count_j + 1; %pass the first calculation into the first row of the storage variables
end
%% Defining the function for the blasius solution
% define initial conditions
function df = Blas_Sc(eta,f,Sc)
df = [f(2); f(3); -0.5*f(1)*f(3);f(5);-0.5*Sc*f(1)*f(5)];
end
