Parametric Continuation Method

A numerical technique used to solve nonlinear equations or systems of equations by gradually varying a parameter.

% Full MATLAB Code to solve F⁽⁴⁾, θ″, and φ″ equations
clear; clc;
% Parameters (example values; update as needed)
k = 1; eta_max = 10; N = 500;
Pr = 0.72; Ra = 1; Ec = 0.1;
A = 1; ab = 0.5;
Br = 0.01; M = 1; S = 1;
Nb = 0.1; Nt = 0.1; Sc = 1; K = 1;

% Mesh
h = eta_max / (N - 1);
eta = linspace(0, eta_max, N);

% Initialize variables
F = zeros(1, N); F1 = zeros(1, N); F2 = zeros(1, N); F3 = zeros(1, N); F4 = zeros(1, N);
theta = ones(1, N); theta1 = zeros(1, N); theta2 = zeros(1, N);
phi = ones(1, N); phi1 = zeros(1, N); phi2 = zeros(1, N);

% Initial conditions
F2(1) = 1;
for i = 2:N
    eta_k = eta(i);
    % F⁽⁴⁾
    F4(i) = (-2*ab*(F2(i-1)^2 + F(i-1)*F3(i-1)) + 2*K*(F3(i-1)*F2(i-1) - F2(i-1)^2/(K+eta_k)) ...
             - A*K*F2(i-1)/(K+eta_k) + K*F(i-1)*(F2(i-1)^2 - F1(i-1)*F3(i-1))/(K+eta_k)) ...
            / (2*ab*(K + eta_k)^2);
    F3(i) = F3(i-1) + h*F4(i-1);
    F2(i) = F2(i-1) + h*F3(i-1);
    F1(i) = F1(i-1) + h*F2(i-1);
    F(i)  = F(i-1)  + h*F1(i-1);

    % θ″
    dF = (F3(i) - F2(i))/(k + eta_k);
    theta2(i) = -theta1(i-1)/(k + eta_k) ...
        - (Ec*Pr/(1+Ra))*(A*dF^2 + ab*Pr*dF^4/(3*(1+Ra))) ...
+ (k*Pr*F(i)*theta1(i-1))/((1+Ra)*(k + eta_k)) ...
        - (Br*Pr*M*F(i)^2)/(1+Ra) - (Pr*S*theta(i-1))/(1+Ra) ...
        - (Nb*Pr*phi1(i-1)*theta1(i-1))/(1+Ra) ...
        - (Nt*Pr*theta1(i-1)^2)/(1+Ra);
    theta1(i) = theta1(i-1) + h*theta2(i-1);
    theta(i) = theta(i-1) + h*theta1(i-1);

    % φ″
    phi2(i) = (Sc*K/(h+K))*F(i)*phi1(i-1) ...
        - (Nt/Nb)*(theta2(i) + theta1(i)/(k + eta_k)) ...
        - phi1(i)/(k + eta_k);
    phi1(i) = phi1(i-1) + h*phi2(i-1);
    phi(i) = phi(i-1) + h*phi1(i-1);
end
% Plots
subplot(3,1,1); plot(eta, F, 'r'); title('F(\eta)');
subplot(3,1,2); plot(eta, theta, 'b'); title('\theta(\eta)');
subplot(3,1,3); plot(eta, phi, 'g'); title('\phi(\eta)');clear; clc;% Parameters (example values; update as needed)k = 1; eta_max = 10; N = 500;Pr = 0.72; Ra = 1; Ec = 0.1;A = 1; ab = 0.5;Br = 0.01; M = 1; S = 1;Nb = 0.1; Nt = 0.1; Sc = 1; K = 1;% Meshh = eta_max / (N - 1);eta = linspace(0, eta_max, N);% Initialize variablesF = zeros(1, N); F1 = zeros(1, N); F2 = zeros(1, N); F3 = zeros(1, N); F4 = zeros(1, N);theta = ones(1, N); theta1 = zeros(1, N); theta2 = zeros(1, N);phi = ones(1, N); phi1 = zeros(1, N); phi2 = zeros(1, N);% Initial conditionsF2(1) = 1;for i = 2:N  eta_k = eta(i);  % F⁽⁴⁾  F4(i) = (-2*ab*(F2(i-1)^2 + F(i-1)*F3(i-1)) + 2*K*(F3(i-1)*F2(i-1) - F2(i-1)^2/(K+eta_k)) ...           - A*K*F2(i-1)/(K+eta_k) + K*F(i-1)*(F2(i-1)^2 - F1(i-1)*F3(i-1))/(K+eta_k)) ...          / (2*ab*(K + eta_k)^2);  F3(i) = F3(i-1) + h*F4(i-1);  F2(i) = F2(i-1) + h*F3(i-1);  F1(i) = F1(i-1) + h*F2(i-1);  F(i)  = F(i-1)  + h*F1(i-1);  % θ″  dF = (F3(i) - F2(i))/(k + eta_k);  theta2(i) = -theta1(i-1)/(k + eta_k) ...      - (Ec*Pr/(1+Ra))*(A*dF^2 + ab*Pr*dF^4/(3*(1+Ra))) ...+ (k*Pr*F(i)*theta1(i-1))/((1+Ra)*(k + eta_k)) ...      - (Br*Pr*M*F(i)^2)/(1+Ra) - (Pr*S*theta(i-1))/(1+Ra) ...      - (Nb*Pr*phi1(i-1)*theta1(i-1))/(1+Ra) ...      - (Nt*Pr*theta1(i-1)^2)/(1+Ra);  theta1(i) = theta1(i-1) + h*theta2(i-1);  theta(i) = theta(i-1) + h*theta1(i-1);  % φ″  phi2(i) = (Sc*K/(h+K))*F(i)*phi1(i-1) ...      - (Nt/Nb)*(theta2(i) + theta1(i)/(k + eta_k)) ...      - phi1(i)/(k + eta_k);  phi1(i) = phi1(i-1) + h*phi2(i-1);  phi(i) = phi(i-1) + h*phi1(i-1);end% Plotssubplot(3,1,1); plot(eta, F, 'r'); title('F(\eta)');subplot(3,1,2); plot(eta, theta, 'b'); title('\theta(\eta)');subplot(3,1,3); plot(eta, phi, 'g'); title('\phi(\eta)');