Численное решение систем дифференциальных уравнений

clc
close all
clear all
 
syms r h(r) h1
h1(r) = piecewise(r<440, 0, 440<=r<650, 175/2, 650<=r<910, (18750-11*r)/290, 910<=r<=970, 35, r>970,0);
h(r) = 2*h1;
%fplot(h1);
%hold on
%fplot(-h1);
%axis([0 1000 -100 100]);
%title('Форма диска');
n = 3000; %об/мин
w = pi*n/30; %рад/с 
Pm = 99*10^6; %Па
P1 = 15*10^6; %Па
E = 2.1*10^11; %Па
mu = 0.3;
ro = 8*10^3; %кг/м^3
r1 = 440;
r2 = 970;
%F = [-(1+mu)/r (1-mu^2)/(E*h*r);
    E*h/r -(1-mu)/r ];
%g = [0; -ro*w^2*h*r];
GU1 = [0; -P1*h(r1)];
space = [440 970];
[r,y] = ode45(@(r,y) SDU1(r,y,mu,E,h,ro,w), space ,GU1);
 
 
function dydr = SDU1(r,y,mu,E,h,ro,w)
dydr1 =-y(1)*(1+mu)/r+y(2)*(1-mu^2)/(E*h*r)-ro*w^2*h*r;
dydr2 = y(1)*E*h/r-y(2)*(1-mu)/r;
dydr = [dydr1;dydr2];
end

function dydr = SDU1(r,y,mu,E,h,ro,w)
h1(r) = piecewise(r<440, 0, 440<=r<650, 175/2, 650<=r<910, (18750-11*r)/290, 910<=r<=970, 35, r>970,0);
if r < 440
    h1 = 0;
elseif 440<=r && r<650
    h1 = 175/2;
elseif 650<=r && r<910
    % ...
end
h = 2*h1;
dydr1 =-y(1)*(1+mu)/r+y(2)*(1-mu^2)/(E*h*r)-ro*w^2*h*r;
dydr2 = y(1)*E*h/r-y(2)*(1-mu)/r;
dydr = [dydr1;dydr2];
end