MATLAB: Using fsolve over a range of equations.

fsolverangesimultaneous

Hi there,
So I have a big array of results, and I want to solve four simultaneous equations over a range.
The code I have sucessfully solves the equations for one case, what I would like to do is run it again for a range of z=100:1:105. I have tried adding loops in and defineing z globally in a loop but none has worked. What I have so far is
This first section is just some reshaping of my external data not too important
function F=simul(c);
READ1=load(['C:\LBM Results\READ1.dat']);
velprof=load(['C:\LBM Results\u-x10000.dat']);
velprof=reshape(velprof,READ1(1),READ1(2));
y1=70;
y2=71;
y3=72;
y4=73;
It is here that I would like z to increment by one once the equations have been solved
z=100;
x1=velprof(z,y1);
x2=velprof(z,y2);
x3=velprof(z,y3);
x4=velprof(z,y4);
F=[-y1+c(1)+c(2)*x1+c(3)*x1^2+c(4)*x1^3;
  -y2+c(1)+c(2)*x2+c(3)*x2^2+c(4)*x2^3;
  -y3+c(1)+c(2)*x3+c(3)*x3^2+c(4)*x3^3;
  -y4+c(1)+c(2)*x4+c(3)*x4^2+c(4)*x4^3];
end
Then to run it simply, this is in a different document to the function above
x0 = [50;0;0;0];
options=optimset('Display','iter'); 
[c,fval] = fsolve(@simul,x0,options)
I hope I have explained it properly!

Best Answer

Solve the equations once and then plug in the values.
for z = 100:105
  x1=velprof(z,y1);
  x2=velprof(z,y2);
  x3=velprof(z,y3);
  x4=velprof(z,y4);
  c = [ (((x3*y4-x4*y3)*x2^2+(-x3^2*y4+x4^2*y3)*x2+x3*x4*y2*(-x4+x3))*x1^3+((-x3*y4+x4*y3)*x2^3+(-x4^3*y3+x3^3*y4)*x2+y2*x3*x4^3-x3^3*x4*y2)*x1^2+((x3^2*y4-x4^2*y3)*x2^3+(x4^3*y3-x3^3*y4)*x2^2+x3^2*x4^2*y2*(-x4+x3))*x1-x3*x2*x4*y1*(-x4+x3)*(x2-x4)*(x2-x3))/((-x4+x3)*(x2-x4)*(x2-x3)*(x1-x4)*(x1-x3)*(x1-x2)), ...
  (((-y2+y1)*x4^2+(y2-y4)*x1^2-x2^2*(y1-y4))*x3^3+((-y1+y2)*x4^3+(y4-y2)*x1^3+x2^3*(y1-y4))*x3^2+((y3-y2)*x1^2+x2^2*(y1-y3))*x4^3+((y2-y3)*x1^3-x2^3*(y1-y3))*x4^2+x1^2*x2^2*(y3-y4)*(x1-x2))/((-x4+x3)*(x2-x4)*(x2-x3)*(x1-x4)*(x1-x3)*(x1-x2)), ...
  (((y4-y3)*x2+(y2-y4)*x3-x4*(y2-y3))*x1^3+((y3-y4)*x2^3+(y4-y2)*x3^3+x4^3*(y2-y3))*x1+((-y1+y4)*x3+x4*(y1-y3))*x2^3+((y1-y4)*x3^3-x4^3*(y1-y3))*x2-x3*x4*(-x4+x3)*(x3+x4)*(-y2+y1))/((-x4+x3)*(x2-x4)*(x2-x3)*(x1-x4)*(x1-x3)*(x1-x2)), ...
  (((y3-y4)*x2+(y4-y2)*x3+x4*(y2-y3))*x1^2+((y4-y3)*x2^2+(y2-y4)*x3^2-x4^2*(y2-y3))*x1+((y1-y4)*x3-x4*(y1-y3))*x2^2+((-y1+y4)*x3^2+x4^2*(y1-y3))*x2+x3*x4*(-y2+y1)*(-x4+x3))/((-x4+x3)*(x2-x4)*(x2-x3)*(x1-x4)*(x1-x3)*(x1-x2)) ];
  disp(z)
  disp(c)
end
Related Question