Why is my code not working?

I am simply plotting a function (temperature) with respect to time(variable).Pardon the lengthy function but i have double checked that there are no errors. I have defined the functions and the deciding constants. i have tried to plot the function with time but resulting in error. I am attaching the code and the questions along with the reference paper. Pls help me.

Answer 1

[Bora Eryilmaz]

The variables (ad, ag, etc.) that you need in your functions T(t), etc. need to be defined in the functions themselves.

Also, fplot() passes a vector of t values to the function T(t), but that function is not written to be handle vector values. You would need to use elementwise multiplications (e.g., t.*t instead of t*t), otherwise you end up with unintended matrices, etc. It is easier to use a loop to compute the values one at a time.

lambdad=56.72; %thermal conductivity of disc
lambdag=1.212; %thermal conductivity of pad
ad=0.00001609; %thermal diffusivity of disc
ag=0.0000003188; %thermal diffusivity of pad
dd=0.016; %thickness of disc
dg=0.01; %thickness of pad
tf=4.72; %braking time
sigma=(lambdag/lambdad)*sqrt(ad/ag);
A=(1-sigma)/(1+sigma);
B=(1-sigma)/((1+sigma)*(1+sigma));
C=((1-sigma)*(1-sigma))/((1+sigma)*(1+sigma)*(1+sigma));

t=0:tf/100:tf+3;
y = zeros(size(t));
for i = 1:numel(t)
  y(i) = T(t(i));
end
plot(t,y)

function Temp=T(t)
lambdad=56.72; %thermal conductivity of disc
lambdag=1.212; %thermal conductivity of pad
ad=0.00001609; %thermal diffusivity of disc
ag=0.0000003188; %thermal diffusivity of pad
dd=0.016; %thickness of disc
dg=0.01; %thickness of pad
tf=4.72; %braking time
sigma=(lambdag/lambdad)*sqrt(ad/ag);
A=(1-sigma)/(1+sigma);
B=(1-sigma)/((1+sigma)*(1+sigma));
C=((1-sigma)*(1-sigma))/((1+sigma)*(1+sigma)*(1+sigma));

Temp=((33040000*sqrt(ad))/(lambdad*(1+sigma)))*(((2*sqrt(t))/sqrt(pi))*(1-((2*t)/(3*tf)))+((4*sqrt(t))/(1+sigma))*(ierfc(2*kd(t))-ierfc(kg(t))-(4*(t/tf)*(inerfc(3,2*kd(t))-(sigma*inerfc(3,2*kg(t))))))+2*(1-(2*A*A))*sqrt(t)*(ierfc(2*(kg(t)+kd(t)))-(4*(t/tf)*(inerfc(3,2*(kg(t)+kd(t)))))+(4*B*sqrt(t))*(ierfc(4*kd(t))-(sigma*ierfc(4*kg(t)))-(4*(t/tf)*(inerfc(3,(4*kd(t)))-(sigma*inerfc(3,(4*kg(t)))))))+2*(1+(3*A)-(A*A))*sqrt(t)*(ierfc(2*((2*kd(t))+kg(t)))-(4*(t/tf)*inerfc(3,2*(2*kd(t)+kg(t)))))+2*(1-(3*A)-(A*A))*sqrt(t)*(ierfc(2*((2*kg(t))+kd(t)))-(4*(t/tf)*inerfc(3,2*(2*kg(t)+kd(t)))))+4*C*sqrt(t)*(ierfc(6*kd(t))-(4*(t/tf)*inerfc(3,6*kg(t))))));
end

function Inerfc=inerfc(n,x)
if n==2
        Inerfc=(erfc(x)-(2*x*(1/sqrt(pi))*(exp(-x*x)-(x*erfc(x)))))/(2*n);
else 
        Inerfc=(((1/sqrt(pi))*(exp(-x*x)-(x*erfc(x))))-(2*x*((erfc(x)-(2*x*(1/sqrt(pi))*(exp(-x*x)-(x*erfc(x)))))/(2*n))))/(2*n);
end
end

function Ierfc=ierfc(x)
Ierfc=(1/sqrt(pi))*(exp(-x*x)-(x*erfc(x)));
end

function Kg= kg(t)
ad=0.00001609; %thermal diffusivity of disc
dd=0.016; %thickness of disc

Kg=dd/sqrt(16*ad*t);
end

function Kd= kd(t)
ag=0.0000003188; %thermal diffusivity of pad
dg=0.01; %thickness of pad

Kd=dg/sqrt(16*ag*t);
end