clear all
d=0.1;TR=10;R=5;
Delta=0.2;
r=0:Delta:R;
n=length(r)-1; %n is the number of intervals
T=zeros(1,n+1);
%Note that I am only solving for the inner variable!
Tinner = fsolve(@(x) HeatedDiskExample(x,r(2:n),R,d,TR,Delta),10*ones(1,n-1));
%Putting the boundary values back in the solution vector.
T(1)=Tinner(1);
T(2:n)=Tinner(1:end);
T(n+1)=TR;
plot(r,T,’ko-‘,’linewidth’,2,’Markersize’,10,’Markerfacecolor’,[0 0 1])
xlabel(‘r’,’FontSize’,14)
ylabel(‘T’,’FontSize’,14)
function f = HeatedDiskExample(x,r,R,d,TR,Delta)
f(1) = (x(2)-2*x(1)+x(1))/Delta^2+(1./(2*Delta*r(1)))*(x(2)-x(1))-d*(x(1).^4-TR.^4)+(r(1)-R).^2;
for i=2:length(r)-1
f(i) = (x(i+1)-2*x(i)+x(i-1))/Delta^2+(1./(2*Delta*r(i)))*(x(i+1)-x(i-1))-d*(x(i).^4-TR.^4)+(r(i)-R).^2;
end
k=length(r);
f(k) = (TR-2*x(k)+x(k-1))/Delta^2+(1./(2*Delta*r(k)))*(TR-x(k-1))-d*(x(k).^4-TR.^4)+(r(k)-R).^2;
end