function walzprozess
clc; 
disp(['Walzprozess',date])
global a b C ha hb Jw Jz Jges mW mZ mges k r1 r2 s0 tspan ys rho B Z
% Konstanten
a=B/2 + Z/2; % Bandmitte - Zapfenmitte [m]
b=0.8; % Bandbreite [m]
C=1850000000; % Geruestmodul [N/m]
Jw=1/4*mW*r1^2+1/12*mW*B^2
Jz=1/4*mZ*r2^2+2*mZ*((B)/2+(Z)/2)^2
Jges=Jw+2*Jz+2*mZ*((B)/2+(Z)/2)^2
mW=pi*r1^2*B*rho %Ballenmasse [kg]
mZ=pi*r2^2*Z*rho %Zapfenmasse [kg]
mges=mW+2*mZ
k=190000000; %Umformwiderstand [N/m^2]
r1=0.5; %Ballenradius [m]
r2=r1/1.5; %Lagerzapfenradius [m]
s0=0.03; %Voreinstellung [m]
rho=7865; %Dichte von Stahl [kg/m^3]
B=b*2; %Ballenbreite
Z=b/1.5; %Zapfenbreite
%Eingangsparameter
tspan=[0;3]
ys=[0;0;0;0];
%DGL Lösen
options=odeset('OutputFcn',@odeplot,'OutputSel',[1 3]);
[t,x]=ode23(@ma,tspan,ys);
%Ausgabe der Ergebnisse
figure, nor=1e+3;
subplot(2,1,1),plot(t,x(:,1)*nor),xlabel('t in [s]'),ylabel('Auffederung in [mm]')
subplot(2,1,2),plot(t,x(:,3)),xlabel('t in [s]'),ylabel('Winkel in [rad]')

function [dgl]=ma(t,x)
global a b C ha hb Jw Jz Jges mW mZ mges k r1 r2 s0 tspan ys rho B Z
%Berechnung der Voreinstellung
ha=0.010*(1+sin(100*t)/5); %Bandhöhe links [m]
hb=0.018*(1+cos(100*t)/5); %Bandhöhe rechts [m]
%DGL
dgl=[x(2);...
    (1/2*sqrt(2)*b*k*sqrt(2)*r1*(ha-x(1)-s0-b*sin(x(3))))-2*C*x(1)+(1/2)*sqrt(2)*b*k*sqrt(2)*r1*(hb-x(1)-s0+b*sin(x(3)))/mges;....
    x(4);...
    (((1/4)*sqrt(2)*b^2*k*sqrt(2)*r1*(ha-x(1)-s0-b*sin(x(3))))-2*C*a*x(1)-(1/4)*sqrt(2)*b^2*k*sqrt(2)*r1*(hb-x(1)-s0+b*sin(x(3))))/Jges];