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

[AKennedy]

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?

[T K]

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);
    a=1;a=1;b=0.1;v=1;u=1;
    M=3;    gamma=pi/3;
Nr=0.1;
     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);
    
    % 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);
        % 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

[T K]

when run this code,the message show .

please help me

[Torsten]

It may happen that "bvp4c" changes the number of discretizations points so that it won't return the solution in 100 points as at the start, but in more or less.

Use

sol = bvp4c(@projfun, @projbc, solinit, options);
Z(i,:) = deval(sol,m,6)

to interpolate the solution to the values m of your choice.

Take care not to come into conflict with the parameter "m" used here:

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))));

[T K]

the m in the relation kh replaced by s s=5.7, Kh=.....