clear

 // Electric parameters
 Ua = 6;     // [V] Armature voltage
 Ra = 21.2;  // [Ohm] Armature resistance at 22°C
 R = Ra;
 La = 217e-6; // [H] Armature inductance
 kE = 4.1157e-3; // [V/rad/s] Back-EMF const.
 // Mechanical parameters
 kM = 4.12e-3; // [Nm/A] Torque const.
 J = 5.2e-9; // [kgm2] Rotor inertia
 Kf = 2.414e-8; // [Nms] Viscous friction
 Tm = 1.13e-3; // [Nm] Stall torque
 Tl = 0; // [Nm] Load torque
 // Thermal parameters
 Rth1 = 20; // [K/W] Thermal resistance
 Rth2 = 48; // [K/W] Thermal resistance
 alpha = 0.004; // [1/K] Temperature coefficient
 Tamb = 22; // [°C] Ambient temperature

 function xdot=DCmotor(t, x)

 w = x(1);
 ia = x(2);

 ui = kE * w;
 Te = kM * ia;
 Ploss = ia^2 * R;
 dTemp = Ploss * (Rth1 + Rth2);
 R = Ra * (1 + alpha*dTemp);
 Twind = dTemp + Tamb;

 xdot(1) = (Te - Tl - Kf*w) / J;
 xdot(2) = (Ua - ui - R*ia) / La;

 endfunction

 taum = 6.5e-3; // [s] Mechanical time const.
 t = 0 : taum/100 : 10*taum;
 t0 = 0;
 x0 = [0; 0];

 x = ode(x0,t0,t,DCmotor);

 figure(1)
 plot(t,x(1,:)*60/2/%pi, "linewidth", 3);
 xgrid(5, 1, 7);
 xlabel("time [s]", "fontsize", 3, "color", "blue");
 ylabel("Speed [rpm]", "fontsize", 3, "color", "blue");
 figure(2)
 plot(t,x(2,:), "linewidth", 3);
 xgrid(5, 1, 7);
 xlabel("time [s]", "fontsize", 3, "color", "blue");
 ylabel("Armature current [A]", "fontsize", 3, "color", "blue");