Solving coupled BVP and ODE without loop

I am looping to solve a second order boundary value problem (1D concentration profile) numerically. This BVP is coupled with an ODE which can be easily solved with the values of the concentration profile solution. There is a 2D mesh required across time and radial position.

The process is as follows:

1. I use bvp4c to solve the boundary value problem at t=0 with the initial condition of the ODE, to find a concentration profile across the radial position.
2. The ODE can then be solved, to find the next time-stepped value, which is required to find the next time-stepped solution for the BVP concentration profile.
3. The time-stepped value from the ODE solution is substituted into the bvp4c solver to find the subsequent profile at the next time step.
4. This is looped through all values of time and space in the mesh.

Obviously this is a time-expensive method, but so far I have not been able to pass a vector through my BVP solver to avoid at least one of the loops. When I attempt to vectorize the BVP function in the bvp4c solver I get errors such as "Dimensions of arrays being concatenated are not consistent." or "The derivative function ODEFUN should return a column vector of length 2".

Any help in reducing processing time would be greatly appreciated!

Here is my current code:

%% Solving 
clear

N1=11; %number of time steps
N2=11; %number of radius steps
N3=10001; %number of mesh grid points in BVP solver
guess = [1; 0];
xmesh = linspace(0,1,10001);
solinit = bvpinit(xmesh,guess);
psimatrix=zeros(N1,N2);
epsmatrix=zeros(N1,N2);
dpsidetamatrix=zeros(N1,N2);


 for j=1:N2
    for i=1:(N1-1)
       
          
epsmatrix(1,:)=1;
eps=epsmatrix(i,j); 

if epsmatrix(i,j)>0
    
f = @(x,y)emdenode(x,y,eps);
S = [0 0; 0 -2]; 
options = bvpset('SingularTerm',S);
sol = bvp4c(f, @emdenbc, solinit, options); %solve
x = linspace(0,1,N3);
dpsidetamatrix(i,j)=sol.yp(1,1+(j-1)*(N3-1)/(N2-1));
psimatrix(i,j)=sol.y(1,1+(j-1)*(N3-1)/(N2-1)); %store solution of psi across radial position in psi matrix
eta=linspace(0,1,N2);
tau=linspace(0,2,N1);
psi=psimatrix(i,j);
eps=1-psi*tau;
epsmatrix(i+1,j)=eps(i+1);
else
    psimatrix(i,j)=1;
epsmatrix(i+1,j)=0;
epsmatrix(i,j)=0;
end
    

       
    end
    
    eps=epsmatrix(N1,j); 
f = @(x,y)emdenode(x,y,eps);
     S = [0 0; 0 -2]; 
options = bvpset('SingularTerm',S);
sol = bvp4c(f, @emdenbc, solinit, options); %solve
x = linspace(0,1,N3);
dpsidetamatrix(N1,j)=sol.yp(1,1+(j-1)*(N3-1)/(N2-1));
psimatrix(N1,j)=sol.y(1,1+(j-1)*(N3-1)/(N2-1)); %store solution of psi across radial position in psi matrix
eta=linspace(0,1,N2);

 end
 

%% Functions



function res = emdenbc(ya,yb)
res = [ya(2)
       yb(1) - 1];
end


function dydx = emdenode(x,y,eps)
sigmahat=0.9977; %sample constants
sigma=sigmahat*sqrt(2*3*3);
%eps=1E-6;

dydx = [y(2) 
       sigma^2*eps.^2*y(1)];

end