How to create a surface for an ODES system that have a variable parameter

Hi everybody.

I like to make a surface for the following system of odes.

the code of the function is this

function [H, h, ZH, Zh, Z] = fun123(muH)
    clear all;
    close all;
    clc;
%
    tic;
    rand('state',sum(100*clock)); % seed
%

    numreps= 1; % Number of iterations
    for j=1:numreps
        options = odeset('RelTol',1e-6,'Stats','on');  
        %Parameters values 
        aH = 800; 
        ah = 0.16; 
        %
        bH =  0.73;
        bh = 0.73;
        %
        alphaH = 1.62;                
        alphah = 1.62;  
        %
        lambdaH = 6.6*10^6; 
        lambdah = 6.6*10^6; 
        %
        rH = 124*52; 
        rh = 124*52; 
        %
        muH = 117*52; %this is the parameter that varying
        muh = 117*52;
        %
        muZ = 45*52; 
        %
        betaH = 6*10^-9;                   
        betah = 6*10^-9; 
        %
        phiH = 1*10^4;                    
        phih = 1*10^4; 
        %
        varrhoH1 = 10^-2;
        varrhoH = varrhoH1*lambdaH; 
        %
        varrhoh1 = 10^-2;                                       
        varrhoh = varrhoh1*lambdah; 
        %
        KH = 1*10^5; % adults host
        Kh = 1*10^5; % tadpoles host
        KZH = 5.6*10^6; % zoospores tadpoles
        KZh = 5.6*10^6; % zoospores tadpoles
    
    
        %x(1) = H; x(2) = h; x(3) = ZH; 
        %x(4) = Zh; x(5) = Z;
     
        G = @(t, x, ah, aH, bh, bH, muh, muH, alphah, ...
        alphaH, lambdah, lambdaH, rh, rH, phih, phiH, ...
        varrhoh, varrhoH, muZ, betaH, betah, Kh, KH, KZH, KZh) ...
        [ah *  x(2) * exp(-Kh * x(2)) - bH * x(1) - alphaH * x(3); ...
        aH * x(1) - (bh + ah * exp(-Kh * x(2))) * x(2) - alphah * x(4); ...
        rH * x(3) * exp(-KZH * x(3)) + lambdaH * x(3) * (x(1)./(varrhoH + x(1))) + x(5) * betaH * (x(1)./(x(1) + x(2))) - x(3) * (bH + muH) - alphaH * x(1) * (x(3)./x(1) + (x(3)./x(1))^2 * ((phiH + 1)./phiH)); ...
        rh * x(4) * exp(-KZh * x(4)) + lambdah * x(4) * (x(2)./(varrhoh + x(2))) + x(5) * betah * (x(2)./(x(1) + x(2))) - x(4) * (bh + ah * exp(-Kh * x(2)) + muh) - alphah * x(2) * (x(4)./x(2) + (x(4)./x(2))^2 * ((phih + 1)./phih)); ...
        lambdaH * x(3) + lambdah * x(4) - x(5) * (muZ + betaH * (x(1)./(x(1) + x(2))) + betah * (x(2)./(x(1) + x(2))))-(lambdaH * x(3) * (x(1)./(varrhoH + x(1))) + lambdah * x(4) * (x(2)./(varrhoh + x(2))))];
        tspan = [0:0.00001:80];
        x0 = [100 100 10 10 500];
        [t,xa] = ode45(@(t,x) G(t, x, ah, aH, bh, bH, muh, muH, ...
        alphah, alphaH, lambdah, lambdaH, rh, rH, phih, phiH, ...
        varrhoh, varrhoH, muZ, betaH, betah, Kh, KH, KZH, KZh), ...
        tspan, x0, options);
        H = xa(:,1)';
        h = xa(:,2)';
        Zh = xa(:,3)';
        ZH = xa(:,4)';
        Z = xa(:,5)';
    end
end

I will like to make a surface that describe the behavior of each solution (H(t), h(t), ZH(t), Zh(t), Z(t)) when a parameter varying, for example \mu_{H}, let say [1000: 500: 6000]. The surface that I want is similar to this figure

Evaluating µH in the range of values given between [1000: 500: 6000]. The code I am implementing to generate the above surface is as follows:

muH = (1*10^3:500:6*10^3);
HMC = [];
ZHMC = [];
fori = 1:length(muH)
    [H, ZH] = fun123(muH(i)); %To evaluate each mu values in the function
    HMC = [HMC; H];
    ZHMC = [ZHMC; ZH];
end

However, I think the result is not correct, since that the \mu_ {H} value does not evaluated correctly within of function.

For example, if I manually evaluate the odes system for a set of values of µ the behavior in H (t) at each value is as follows

If \mu=1000, the graphic of H(t) in 2-D is this

If \mu=2000, the graphic of H(t) in 2-D is this

If \mu=3000, the graphic of H(t) in 2-D is this

In this sense, I hope someone else can I help me for make this surface in 3-D for the range values of \mu

Answer 1

[Ameer Hamza]

Inside the function definition, You are again setting it to a constant value

muH = 117*52; %this is the parameter that varying
muh = 117*52;

No matter what you pass as muH input, it gets overwritten and you same constant output. You should remove this line from the function definition.

[Alvaro Sepulveda]

But, I have this error

[Ameer Hamza]

I think you are directly trying to run your function without passing the value of muH. Are you running this script?

muH = (1*10^3:500:6*10^3);
HMC = [];
ZHMC = [];
fori = 1:length(muH)
    [H, ZH] = fun123(muH(i)); %To evaluate each mu values in the function
    HMC = [HMC; H];
    ZHMC = [ZHMC; ZH];
end

[Alvaro Sepulveda]

Yes, with this script I want to evaluate the values of muH and obtain a curve of H (t) for each value of muH and generate the surface above

[Ameer Hamza]

Oh! I missed it initially. Remove these lines from your function definition

clear all;
close all;

[Alvaro Sepulveda]

Hi Ameer Hamza, thank you. I can do it! I removed the last lines above.

[Ameer Hamza]

Does the plot looks correct now?

[Alvaro Sepulveda]

Hi Ameer Hamza,

Yeah!

Look how it is for example ZH(t), Thanks for your suggestions!

[Ameer Hamza]

You can run

shading interp

to remove the black lines.

[Alvaro Sepulveda]

Ummm.. ok, I need to make other graphics and if I have anything doubt I will write a new post, very thanks Ameer Hamza

regards