Finitadiffrential[numeriska metoder-MATLAB]

clc,clf,clear all
%vår info
u0=450; r0=1;R=2;K=1;uom=20;
y=1:2:20;
x=2*y;  %används för att fördubbla våra intervall
g=0;  %fortsätt läs koden!
for i=1:numel(x)
    
    N=25*x(i);  %fördubblar delintervallena
    n=N;
    h=(R-r0)/N;  %steglängd
    r=r0:h:R;
    a=((8*(r.^3)-27*(r.^2)+24*r-6)./(2*(r.^3)-9*(r.^2)+12*r-6));
    sup=(2.*r+h.*a)';  %super diagonal
    dia=(-4.*r)';       %huvud diagonal
    sub=(2.*r-h.*a)';  %sub diagonal
    
    %Vår vänsterled------
    A=spdiags([sup dia sub],-1:1,N,N);% bildar en gles matris
    A(end,end-1)=-(2*r(end)+(h^2)*a(end)*K);  % byter ut den första elementen 
    % med vårt handberäknade värde
    A(end,end)=r(end);  % byter ut den sista elementen
    A(end,end-2)=r(end);
    %Vår högerled
    
    B=zeros(N,1);  % högerled
    B(1,1)=-u0*(2*r(1)-h*a(1));  % byter ut elementen med 
    % vårt handberäknade värde
    B(end,1)=-20*a(end)*(h^2)*K;  % byter ut elementen 
    %Beräknar 
    
    U=A\B;
    u=[u0 U'];
    U=u(end);
    
    %visar resultatet
    disp(['temperaturen vid r=2 med ',num2str(N),' delintervall är: ',num2str(u(end))])
    %Ritar vår graf
    g=g+1;    %används för att rita olika diagram
   
    plot(r,u),grid on,hold on
    xlabel('Radie')
    ylabel('Temperatur')
    title(['Diagram med ',num2str(N),' delintervall' ])
end