I was wondering how coupled ODEs could be solved with Matlab.. I'm having trouble passing the coupling solution vectors between solvers since the length of the vectors isn't constant during the iteration e.g. sol1 = bvpinit(linspace(0.,h,200),[0.1 0]); sol2 = bvpinit(linspace(0.,h,200),@mat4init); f=sol1.y g=sol2.y sol1 = bvp4c(@fODE,@fBC,sol1,options1,g,Re,K); sol2 = bvp4c(@gODE,@gBC,sol2,options1,f,Re,K); function dydx = fODE(x,y,g,Re,K) dydx = [ y(2,:) Re*(y(1,:).^2-g.*2+K^2)]; function dgdx = gODE(x,y,g,Re,K) dgdx = [ g(2,:) 2*Re*(y(1,:).*g(1,]; this doesn't work bacause the dimensions of y and g dont match during iteration. Any ideas how to overcome this?