This example works on FinSim 10_01_31 and subsequent versions.

This example models the force of a tennis ball hitting a wall, based on ideas presented by S J Haake, M J Carre, R Kirk, and S R Goodwill
at the University of Sheffield, Dept. of Mechanical Engineering..
The main model is that of a spring-mass system governed by the formula:
m*y(2) + c*y(1) + k*y = 0, with y = 0 and y(1) = speed of ball hitting the wall.
There are three forces that push against the wall:

1) the spring force, spring_force is due to the ball being compressed,
with spring_force = k*y.

2) the damper force, damper_force is due to the viscuosity of the ball
with damper_force = c*y(1)

3) the flux force, flux_force is due to the portion of the ball that loses it's speed during a small time slice. This flux_force is m0 * speed / delta_time, where m0 = m*s1/s2, where s1 is approx. the surface of a cillinder having as base the flattened portion of the ball and as hight the current speed y0 * delta_time and s2 is the surface of the ball, which makes

flux_force = y0*y0*m* $VpSqrt(2*r*x0-x0*x0)/(2*r*r), where y0 is the speed at a particular time, m is the total mass of the ball, r is the radious of the ball, and x0 is the position of the center of the ball with respect to its original position when the impact began.

The algorithm solves the problem for speeds of 10m/s and 30m/s in a loop. In an inner loop the algorithm adjusts the non-linear coefficients and solves the equation for a number of steps with constant coeficients. It then performs again the same steps for the new value of the loop variable (named slice), which automatically adjusts which portion of the data is provided to the differential eq. solver via the "view..as" mechanism. Note that by using the "view..as" mechanism no data transfer is needed during the computation, the data being computed "in place".

The results are displayed via the $PrintM task as well as via the $VpPtPlot task. The latter produces the results shown here

module top;
parameter real alpha = 1.65;
parameter real r_total_time = 0.005;/* seconds */
parameter nr_pts_per_ct_coef = 10;
parameter nr_slices_per_sec = 10000;
parameter nr_pts_per_sec = nr_pts_per_ct_coef * nr_slices_per_sec;
parameter integer Size = r_total_time*nr_pts_per_sec;
parameter real r_size = Size;
parameter real h = 2*r_total_time/r_size;/* set double of sampling period */
parameter real m = 0.057;/* kg */
parameter real A = 16000000;/* N/m**2 */
parameter real k0= 21000; /*N/m*/
parameter real B = 3500; /* Ns/m */
parameter real r = 0.032; /*m (radius of ball)*/
parameter order = 2;
parameter nrEq = 1;
localparam integer nr_slices = nr_slices_per_sec * r_total_time;
real Fe[0:0][0 : Size], x[0:0][0 : Size], y[0:0][0 : Size];
reg [0 : 3199] ressymb[0 : 0];
real coef[0 : 0][0 : 2];
view real x_ct [0:0][0:nr_pts_per_ct_coef]
as x[0][$I2+slice*nr_pts_per_ct_coef];
view real Fe_ct [0:0][0:nr_pts_per_ct_coef]
as Fe[0][$I2+slice*nr_pts_per_ct_coef];
view real y_ct [0:0][0:nr_pts_per_ct_coef]
as y[0][$I2+slice*nr_pts_per_ct_coef];
real ar_f_spring [0:nr_slices];
real ar_f_damper [0:nr_slices];
real ar_f_flux [0:nr_slices];
real ar_f_total [0:2][0:nr_slices];
integer i, j, isZero, total_time;
integer slice = 0;
real k, c, x0, y0, v0;
initial begin
#1;
coef[0][0] = m;
$InitM(Fe, 0.0);
for (j = 0; j < 2; j = j + 1)
begin
/* find total force for two speeds: v0 = 10m/s and v0 = 30m/s */
v0 = (j == 0) ? 10 : 30;
x[0][0] = 0;
y[0][0] = v0;
isZero = 0;/* kludge for nulifying results for when the model does not apply*/
for (slice = 0; slice <= nr_slices; slice = slice + 1)
begin
x0 = (x_ct[0][0] > 0) ? x_ct[0][0] : -x_ct[0][0];
y0 = y_ct[0][0];
c = 4.0*B*x0*(2.0*r-x0);
k = (slice < 2) ? 120000 : (k0 + A*(x0**alpha));
coef[0][1] = c;
coef[0][2] = k;
ar_f_spring[slice]= x0*k;
ar_f_damper[slice] = y0*c;
if ((y0 > 0) && (x0 > 0)) begin
ar_f_flux[slice] = y0*y0*m*
$VpSqrt(2*r*x0-x0*x0)/(2*r*r);
end
else ar_f_flux[slice] = 0;
ar_f_total[j][slice] = ar_f_spring[slice] +
ar_f_damper[slice] +
ar_f_flux[slice];
ar_f_total[j][slice] = (isZero) ? 0: ((ar_f_total[j][slice]>0)?ar_f_total[j][slice]:0);
if ((slice != 0) && (ar_f_total[j][slice] < 15
/* after this value the spring model no longer applies*/ )) begin
isZero = 1;
end
/* call dif eq solver. The array x will contain the solution and the
array y will contain the first derivative of x.
Note that the initial conditions have been already placed in x[0]
and y[0].
*/
$VpLODE(order, nrEq, h, nr_pts_per_ct_coef+1,
x_ct, coef, Fe_ct, y_ct, ressymb);
end
ar_f_total[j][nr_slices] = 0.0;
end
$PrintM(ar_f_total,"%e");
total_time = r_total_time*1000;
$Flot("ball_force.html", 2, h/2,
"Tennisball_(0.057kg,0.032m)_Force_pushing_the_wall)",
"Time (ms)", "Total Force(N)", 0, nr_slices, ar_f_total,
"10m/s", "30m/s");
end
endmodule