how do i change an oscillatory plot

I have a plot of an analytical and numerical solution of a 2nd order differential equation which should match up on the same path. I'm getting the analytical as the correct plot but, the numerical solution plots out as oscillatory and not on the same path.

Here's the code i used and the output graph for the equation, The problem is given by:

wc * dT/dz - (piD^2/4*k) * d^2T/dz^2 - qs + (pi*Dh)(T-Ta)=0

the analytical solution is for when h = 0

function [ z, T ] = heattransfer( a, b, c , g )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Inputs                                                                 %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% problem parameters: ====================================================

% coefficients:
a = 40.0;      % a = wc
b = 3.65;      % b = (pi*D^2/4) * k;
c = 10000.0;   % c = qs
g = 239.3;     % g = (pi*D*h)*(T-Ta)

% domain:
L = 1.0;

% boundary conditions:
T0 = 30;
TL = 200;   %TL = Ttarget

% solver parameters: =====================================================

% solution domain:
N = 200;               % number of nodes
z = linspace( 0 , L, N ); % set z equally spaced over the z domain

% initial guess: (constant values)
  
Tinit = T0;             %T(z) = T0;
dTdzinit = 0;             %dT/dz(z) = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Solution                                                               %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% solution initialization:
solinit = bvpinit( z, [ Tinit; dTdzinit ] );

% solution process - bvp4c solver:
% - the "@" notation is used in function handles; the expression:
%   "@( var1, var2, ... )fun( var1, var2, ..., vari, ... varn )"
%   allows using variables var1, var2, ... in calls to the function "fun" 
%   while the values of variables vari ... varn specified prior to the
%   call to bvp4c
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)


tic  % start clock for computing time

sol = bvp4c( @( z, y )odefun( z, y, a, b, c, g ), ...
             @( ya, yb )bcfun( ya, yb, T0, TL ), ...
             solinit );
         
toc  % end clock for computing time


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Post-Processing                                                        %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% retrieve solution:

y = deval( sol, z );
T = y( 1, : );


% analytical solution for h = 0:

Tanalytic = T0 + (TL - T0) * ( ( exp( a * z/b ) - 1 )/( exp( a * L/b ) - 1 ) ) ...
    + (c * L/b ) * ( z/L - ( exp( a * z/b ) - 1 )/( exp( a * L/b ) - 1 ) );

% plots:
figure;
plot( z, T, '-b', z, Tanalytic, '+r')
grid;
box on;
legend( ' T ', ' t_{analytic} (for h = 0) ' )
title(' Boundary Value Problem ' )
xlabel( ' z ' )
ylabel( ' T(z) ' )

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Function describing the BVP: dy/dx = F( x, y )                         %
%                                                                        %
% For this problem: y = [ y1; y2 ], y1 = T, y2 = dT/dz                   %
%                                                                        %
%          wc * dT/dz - (piD^2/4*k) * d^2T/dz^2 - qs + (pi*Dh)(T-Ta)=0            %
%                                                                        % 
% => F  = [ F1; F2 ]:                                                    %
%    F1 = dy1/dz = y2 = dT/dz                                            %
%    F2 = dy2/dz = d^2T/dz^2 = ( wc * y2 - qs * y1 - (piDh)(T-Ta) )/(piD^2/4k)        %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dydz = odefun( z, y, a, b, c, g )

%piD24k = pi * D^2/ (4*k)

dydz = [ y(2);
         ( a * y(2) - c * y(1) - g ) / b ];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Function evaluating the boundary conditions: res = F( ya, yb, ... )    %
%                                                                        %
% For this problem: y = [ y1 y2 ], y1 = T, y2 = dT/dz                    %
%                                                                        %
% x = 0: T - T0 = 0 --> y1 - T0 = 0                                      %
% x = L: T - Ttarget = 0 --> y1 - Ttarget = 0                                      %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function res = bcfun( ya, yb, T0, TL )

res = [ ya(1) - T0  ;
        yb(1) - TL ];
    
end

[David Goodmanson]

Hello izuchukwu

Ordinarily the heat conduction equation would go like

wc * dT/dt - (piD^2/4*k) * d^2T/dz^2 - qs + (pi*Dh)(T-Ta)=0
         ^

instead of

wc * dT/dz - (piD^2/4*k) * d^2T/dz^2 - qs + (pi*Dh)(T-Ta)=0
         ^

Is your equation in error, or is there a reason that you were able to substitute d/dz for d/dt ?

[izuchukwu morah]

that's exactly how the original equation from the problem i'm trying to solve was given.