Please help me to run surf plot with bvp4c.

Please help me to run surf plot with bvp4c.The surfce digram consists of (constant Prf [2 :1:6] represents y-axis & vector>> sol.x [0 4] represents x-axis (m = linspace(0,4);) & second solution only sol.y(6,:) represents z-axis).The following is the code for 2D (sol.x [0 4] and only sol.y(6,:)). How to give command for making surf plot in bvp4c.

proj()

rr = 1x4

4 5 6 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The solution was obtained on a mesh of 227 points.
The maximum residual is  8.171e-05. 
There were 7965 calls to the ODE function. 
There were 122 calls to the BC function. 
    0.2422

The solution was obtained on a mesh of 179 points.
The maximum residual is  2.747e-05. 
There were 7696 calls to the ODE function. 
There were 149 calls to the BC function. 
    0.1666

The solution was obtained on a mesh of 188 points.
The maximum residual is  2.582e-05. 
There were 7831 calls to the ODE function. 
There were 149 calls to the BC function. 
    0.1666

The solution was obtained on a mesh of 198 points.
The maximum residual is  3.095e-05. 
There were 7981 calls to the ODE function. 
There were 149 calls to the BC function. 
    0.1935

ans = struct with fields:
    solver: 'bvp4c'
         x: [0 0.0202 0.0404 0.0808 0.1212 0.1818 0.2424 0.3030 0.3636 0.4242 0.4848 0.5253 0.5657 0.6061 0.6465 0.7071 0.7677 0.8081 0.8485 0.8889 0.9293 ... ] (1x198 double)
         y: [9x198 double]
        yp: [9x198 double]
     stats: [1x1 struct]

function sol= proj
clc;clf;clear;
%Relation of base fluid
rhof=1;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;sigf=0.05*10^-8;alfaf=kf/(rhof*cpf);
%FE3O4
ph1=0.01;rho1=5180*10^-3;cp1=670*10^4;k1=9.7*10^5;sig1=0.74*10^-2;

%copper
ph2=0.01;rho2=8933*10^-3;cp2=385*10^4;k2=401*10^5;sig2=5.96*10^-1;

%Relation of hyprid
m=5.7;
kh=kf*((k1+(m-1)*kf-(m-1)*ph1*(kf-k1))/((k1+(m-1)*kf+ph1*(kf-k1))))*((k2+(m-1)*kf-(m-1)*ph2*(kf-k2))/((k2+(m-1)*kf+ph2*(kf-k2))));
muh= muf/((1-ph1)^2.5*(1-ph2)^2.5);
rhoh=rhof*(1-ph2)*((1-ph1)+ph1*(rho1/rhof))+ph2*rho2;
v1 =rhof*cpf*(1-ph2)*((1-ph1)+ph1*((rho1*cp1)/(rho2*cp2)))+ph2*(rho2*cp2);
sigh=sigf+(3*((ph1*sig1+ph2*sig2)-sigf*(ph1+ph2))/(((ph1*sig1+ph2*sig2)/(sigf*(ph1+ph2)))+2-((ph1*sig1+ph2*sig2)/sigf)+(ph1+ph2)));
alfah=kh/v1;
myLegend1 = {};
rr = [4 5 6 7]
for i =1:numel(rr)
    Prf = rr(i);
Nr=0.1;
    gamma=pi/3;
a=1;b=0.1;v=1;u=1;
   M=3;
  Nt=1;Nb=1; sc=0.6;s1=1;s2=1;
  p=-0.5; L=(muf/rhof);L1=L^(p);

   Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
   Lf=rhof*kf;

   y0 = [1,0,1,0,0,1,0,1,0];
options =bvpset('stats','on','RelTol',1e-4);
    m = linspace(0,4);
   solinit = bvpinit(m,y0);
      sol= bvp4c(@projfun,@projbc,solinit,options);
      disp((sol.y(1,20)))
    figure(1)
    plot(sol.x,(sol.y(6,:)))
   % axis([0 4 0 1])
    grid on,hold on
    myLegend1{i}=['Pr = ',num2str(rr(i))];
    

    
    i=i+1;
end
figure(1)
legend(myLegend1)
hold on


    function dy= projfun(~,y)
        dy= zeros(9,1);
%         alignComments
        E  = y(1);
        dE = y(2);
        F  = y(3);
        dF= y(4);
        W = y(5);
        t = y(6);
        dt  = y(7);
        phi = y(8);
        dphi  = y(9);
       
    dy(1) = dE;        
          dy(2) = (rhoh/muh)*((((a*u)/L1^(2)))*E^2+(1/L1)*W*dE+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*E*sin(gamma)*sin(gamma));
         dy(3) = dF;
          dy(4) = (rhoh/muh)*((((b*v)/L1^(2)))*F^2+(1/L1)*W*dF+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*F*sin(gamma)*sin(gamma));

         dy(5) = -(1/L1)*(u*a*E+b*v*F);
        dy(6) = dt;
dy(7) =(Prf*(rhof/muf))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2);
        dy(8)= dphi;
dy(9)=(sc/L^(p+1))*W*dphi-(s1/s2)*(Nt/Nb)*(((Prf*(rhof/muf)))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2));
    
    end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(3)-1; 
    ya(5)-0;
    ya(6)-1;
    ya(8)-1;
    yb(1);
    yb(3);
    yb(6);
    yb(8)];
end

Answer 1

[Abhinaya Kennedy]

To create a surface plot using the bvp4c solution in MATLAB, you need to organize your data into matrices that represent the x, y, and z axes. In your case, the x axis is represented by sol.x, the y axis by Prf, and the z axis by sol.y(6,:). Here's how you can modify your code to produce a surface plot:

1. Initialize Data Storage: Store the results for each value of Prf in a matrix.
2. Loop Over Prf: Solve the boundary value problem for each value of Prf.
3. Create the Surface Plot: Use surf to plot the data.

proj()

Unrecognized function or variable 'dt'.

Error in solution>proj/projfun (line 56)
dt

Error in bvparguments (line 96)
testODE = ode(x1,y1,odeExtras{:});

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

Error in solution>proj (line 34)
sol = bvp4c(@projfun, @projbc, solinit, options);

function sol = proj
    clc; clf; clear;
    % Define constants and parameters
    rhof = 1; kf = 0.613e5; cpf = 4179e4; muf = 1e-3 * 10; sigf = 0.05e-8;
    alfaf = kf / (rhof * cpf);
    ph1 = 0.01; rho1 = 5180e-3; cp1 = 670e4; k1 = 9.7e5; sig1 = 0.74e-2;
    ph2 = 0.01; rho2 = 8933e-3; cp2 = 385e4; k2 = 401e5; sig2 = 5.96e-1;
    m = 5.7;
    kh = kf * ((k1 + (m - 1) * kf - (m - 1) * ph1 * (kf - k1)) / ((k1 + (m - 1) * kf + ph1 * (kf - k1)))) * ...
         ((k2 + (m - 1) * kf - (m - 1) * ph2 * (kf - k2)) / ((k2 + (m - 1) * kf + ph2 * (kf - k2))));
    muh = muf / ((1 - ph1)^2.5 * (1 - ph2)^2.5);
    rhoh = rhof * (1 - ph2) * ((1 - ph1) + ph1 * (rho1 / rhof)) + ph2 * rho2;
    v1 = rhof * cpf * (1 - ph2) * ((1 - ph1) + ph1 * ((rho1 * cp1) / (rho2 * cp2))) + ph2 * (rho2 * cp2);
    sigh = sigf + (3 * ((ph1 * sig1 + ph2 * sig2) - sigf * (ph1 + ph2)) / ...
          (((ph1 * sig1 + ph2 * sig2) / (sigf * (ph1 + ph2))) + 2 - ((ph1 * sig1 + ph2 * sig2) / sigf) + (ph1 + ph2)));
    alfah = kh / v1;

    % Initialize parameters
    rr = [4 5 6 7];
    numPrf = numel(rr);
    m = linspace(0, 4);
    y0 = [1, 0, 1, 0, 0, 1, 0, 1, 0];
    options = bvpset('stats', 'on', 'RelTol', 1e-4);
    
    % Initialize storage for surface plot
    Z = zeros(numPrf, length(m));
    
    % Solve the BVP for each Prf
    for i = 1:numPrf
        Prf = rr(i);
        solinit = bvpinit(m, y0);
        sol = bvp4c(@projfun, @projbc, solinit, options);
        Z(i, :) = sol.y(6, :); % Store the z-axis data
    end
    
    % Create surface plot
    [X, Y] = meshgrid(m, rr);
    figure;
    surf(X, Y, Z);
    xlabel('x');
    ylabel('Prf');
    zlabel('Solution y(6,:)');
    title('Surface Plot of Solution');
    grid on;

    function dy = projfun(~, y)
        dy = zeros(9, 1);
        E = y(1);
        dE = y(2);
        F = y(3);
        dF = y(4);
        W = y(5);
        t = y(6);
        dt
    end
end

[Walter Roberson]

Your function projfun is incomplete: it needs to calculate the current components for dy.

[T K]

proj() 

rr = 1x4

4 5 6 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The solution was obtained on a mesh of 227 points.
The maximum residual is  8.171e-05. 
There were 7965 calls to the ODE function. 
There were 122 calls to the BC function. 
    0.2422

The solution was obtained on a mesh of 179 points.
The maximum residual is  2.747e-05. 
There were 7696 calls to the ODE function. 
There were 149 calls to the BC function. 
    0.1666

The solution was obtained on a mesh of 188 points.
The maximum residual is  2.582e-05. 
There were 7831 calls to the ODE function. 
There were 149 calls to the BC function. 
    0.1666

The solution was obtained on a mesh of 198 points.
The maximum residual is  3.095e-05. 
There were 7981 calls to the ODE function. 
There were 149 calls to the BC function. 
    0.1935

ans = struct with fields:
    solver: 'bvp4c'
         x: [0 0.0202 0.0404 0.0808 0.1212 0.1818 0.2424 0.3030 0.3636 0.4242 0.4848 0.5253 0.5657 0.6061 0.6465 0.7071 0.7677 0.8081 0.8485 0.8889 0.9293 ... ] (1x198 double)
         y: [9x198 double]
        yp: [9x198 double]
     stats: [1x1 struct]

%full code is showen below...
function sol= proj
clc;clf;clear;
%Relation of base fluid
rhof=1;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;sigf=0.05*10^-8;alfaf=kf/(rhof*cpf);
%FE3O4
ph1=0.01;rho1=5180*10^-3;cp1=670*10^4;k1=9.7*10^5;sig1=0.74*10^-2;
%copper
ph2=0.01;rho2=8933*10^-3;cp2=385*10^4;k2=401*10^5;sig2=5.96*10^-1;
%Relation of hyprid
m=5.7;
kh=kf*((k1+(m-1)*kf-(m-1)*ph1*(kf-k1))/((k1+(m-1)*kf+ph1*(kf-k1))))*((k2+(m-1)*kf-(m-1)*ph2*(kf-k2))/((k2+(m-1)*kf+ph2*(kf-k2))));
muh= muf/((1-ph1)^2.5*(1-ph2)^2.5);
rhoh=rhof*(1-ph2)*((1-ph1)+ph1*(rho1/rhof))+ph2*rho2;
v1 =rhof*cpf*(1-ph2)*((1-ph1)+ph1*((rho1*cp1)/(rho2*cp2)))+ph2*(rho2*cp2);
sigh=sigf+(3*((ph1*sig1+ph2*sig2)-sigf*(ph1+ph2))/(((ph1*sig1+ph2*sig2)/(sigf*(ph1+ph2)))+2-((ph1*sig1+ph2*sig2)/sigf)+(ph1+ph2)));
alfah=kh/v1;
myLegend1 = {};
rr = [4 5 6 7]
for i =1:numel(rr)
    Prf = rr(i);
    Nr=0.1;
    gamma=pi/3;
    a=1;b=0.1;v=1;u=1;
    M=3;
    Nt=1;Nb=1; sc=0.6;s1=1;s2=1;
    p=-0.5; L=(muf/rhof);L1=L^(p);
    Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
    Lf=rhof*kf;
    y0 = [1,0,1,0,0,1,0,1,0];
    options =bvpset('stats','on','RelTol',1e-4);
    m = linspace(0,4);
    solinit = bvpinit(m,y0);
    sol= bvp4c(@projfun,@projbc,solinit,options);
    disp((sol.y(1,20)))
    figure(1)
    plot(sol.x,(sol.y(6,:)))
    % axis([0 4 0 1])
    grid on,hold on
    myLegend1{i}=['Pr = ',num2str(rr(i))];
    
    
    i=i+1;
end
figure(1)
legend(myLegend1)
hold on
    function dy= projfun(~,y)
        dy= zeros(9,1);
        % alignComments
        E = y(1);
        dE = y(2);
        F = y(3);
        dF= y(4);
        W = y(5);
        t = y(6);
        dt = y(7);
        phi = y(8);
        dphi = y(9);
        
        dy(1) = dE; 
        dy(2) = (rhoh/muh)*((((a*u)/L1^(2)))*E^2+(1/L1)*W*dE+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*E*sin(gamma)*sin(gamma));
        dy(3) = dF;
        dy(4) = (rhoh/muh)*((((b*v)/L1^(2)))*F^2+(1/L1)*W*dF+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*F*sin(gamma)*sin(gamma));
        dy(5) = -(1/L1)*(u*a*E+b*v*F);
        dy(6) = dt;
        dy(7) =(Prf*(rhof/muf))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2);
        dy(8)= dphi;
        dy(9)=(sc/L^(p+1))*W*dphi-(s1/s2)*(Nt/Nb)*(((Prf*(rhof/muf)))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2));
        
    end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
    ya(3)-1; 
    ya(5)-0;
    ya(6)-1;
    ya(8)-1;
    yb(1);
    yb(3);
    yb(6);
    yb(8)];
end

[Walter Roberson]

Okay... but you do not produce a surface plot?