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.008; %Voreinstellung [m]
rho=7865; %Dichte von Stahl [kg/m^3]
B=b*2; %Ballenbreite
Z=b/1.5; %Zapfenbreite
%Eingangsparameter
tspan=[0;6]
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.02; %Bandhöhe links [m]
hb=0.02; %Bandhöhe rechts [m]
%DGL
dgl=[x(2);...
    ((b*k*sqrt(r1*(ha-s0-2*x(1)+b*sin(x(3))))-2*C*x(1)+b*k*sqrt(r1*(hb-s0-2*x(1)-b*sin(x(3)))))/mges);...
    x(4);...
    (((1/2)*b^2*k*sqrt(r1*(ha-s0-2*x(1)-b*sin(x(3))))-(1/2)*b^2*k*sqrt(r1*(hb-s0-2*x(1)+b*sin(x(3))))-2*C*a^2*sin(x(3)))/Jges)];