Algorithm to find by attempts, with R, the function of a certain type that "better" approximates certain data (the sum of the squares of the differences between the data (input, output) and the images of the inputs is calculated using the function).
    It is necessary to choose a wide initial interval [m1,M1] for the unknown parameter U (or 2 or ... or 6 wide initial intervals [m1,M1], [m2,M2], ..., [m6,M6] for the unknown parameters U, V, W, Z, J, K) which certainly contain the solutions (it is necessary to choose them taking into account the known points and the type of function with which they are to be approximated), and after the command F_RUN1 (or F_RUN2, F_RUN3, ... F_RUN6). Then you can gradually narrow down the intervals.
    You need to choose a number of steps that is higher the more the parameters are (if the parameters are 1,2, or 3, 1e5 is enough; if they are more numerous, it is better to use 1e6 or 2e6 steps).
    Obviously the algorithm must be modified according to the type of function to be searched.
    To use the program first load source("http://macosa.dima.unige.it/r.R") and then the following lines, after manipulating the ones in blue. See examples one, ..., seven.

x = c( ... ); y = c( ... )
# Choose and define the type of function ...
yy1 = function(U,x) ...; Fun = function(x) ... # with U0 instead of U
# or:
yy2 = function(U,V,x) ...; Fun = function(x) ... # with U0,V0 instead of U,V
# or:
yy3 = function(U,V,W,x) ...; Fun = function(x) ...
# or:
yy4 = function(U,V,W,Z,x) ...; Fun = function(x) ...
# or:
yy5 = function(U,V,W,Z,J,x) ...; Fun = function(x) ...
# or:
yy6 = function(U,V,W,Z,J,K,x) ...; Fun = function(x) ... # with U0 ... K0 instead of U ... K
# By F_RUN1 (… F_RUN6), in U0 (V0,W0,Z0) the values that minimize the deviation
# [measured as the sum of the squares of the differences between the ordinates and
# the corresponding values of Fun(x)] are set, and in dif the resulting deviation
# is put, whose square root is printed.
# At the beginning I choose a big deviation, which will be minimized:
dif = 1e300
# I start by choosing intervals where the unknown parameters are for sure:
m1=…; M1=…; F_RUN1(x,y,1e5); cat(U0," d:",NUMBER(dif^0.5),"\n")
# or:
m1=;M1=; m2=;M2=; F_RUN2(x,y,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# or:
m1=;M1=; m2=;M2=; m3=;M3=;  F_RUN3(x,y,1e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# or:
m1=;M1=; m2=;M2=; m3=;M3=; m4=;M4=; F_RUN4(x,y,1e5); cat(U0,V0,W0,Z0," d:",NUMBER(dif^0.5),"\n")
# or:
m1=;M1=; m2=;M2=; m3=;M3=; m4=;M4=; m5=;M5=; F_RUN5(x,y,1e5); cat(U0,V0,W0,Z0,J0," d:",NUMBER(dif^0.5),"\n")
# or:
m1=;M1=; m2=;M2=; m3=;M3=; m4=;M4=; m5=;M5=; m5=;M5=; F_RUN5(x,y,1e5); cat(U0,V0,W0,Z0,J0,K0," d:",NUMBER(dif^0.5),"\n")
# I continue narrowing the intervals; eventually increase the n. of the tests
# ...
# When I want, I draw the graph of the function gradually obtained.
# At the end I draw the graph of the function on which we have stabilized
# If you want more digits use  more( c(U0,V0,...) )

Examples

One
I want to approximate the relationship between two random variables G1 and G2 with a relation of the type G2 = k·G1 + h.
I have to determine k and h. In three experiments I get the pairs (1.16, 18), (3.36, 31), (4.8, 48).

G1 = c(1.6, 3.6, 4.8)
G2 = c(18,  31,  48)
range(G1); range(G2)
# 1.6  4.8    18  48
Plane(0,5, 0,50); POINT(G1,G2,"brown")
yy2 = function(U,V,x) U+V*x; Fun = function(x) U0+V0*x
dif = 1e300
m1=-10; M1=10; m2=0; M2=50; F_RUN2(G1,G2,1e5); cat(U0,V0, " d:",NUMBER(dif^0.5),"\n")
#2.06595 9.074395  d: 4.64685263153777
# let's do the chart to see what we got
graph1(Fun,0,5, "red")
# In this case (linear trend) if I continue further I do not improve significantly:
m1=1; M1=3; m2=8; M2=10; F_RUN2(G1,G2,1e5); cat(U0,V0, " d:",NUMBER(dif^0.5),"\n")
#2.068302 9.079689  d: 4.64670394481248
m1=2; M1=2.08; m2=9; M2=9.1; F_RUN2(G1,G2,1e5); cat(U0,V0, " d:",NUMBER(dif^0.5),"\n")
#2.061353 9.081599  d: 4.64670170566259    2.061 + 9.082*x
# if I want to continue:
m1=2.06; M1=2.065; m2=9.075; M2=9.085; F_RUN2(G1,G2,1e5); cat(U0,V0, " d:",NUMBER(dif^0.5),"\n")
#2.061236 9.081633  d: 4.64670170498558
m1=2.061; M1=2.062; m2=9.081; M2=9.082; F_RUN2(G1,G2,1e5); cat(U0,V0, " d:",NUMBER(dif^0.5),"\n")
#2.061222 9.081633  d: 4.6467017049404     2.0612 + 9.0816*x
graph1(Fun,0,5, "brown")


# Obviously in this case I could have used the regression technique:
regression1(G1,G2)
# 9.082 * x + 2.061
Two
I know that a quantity y is linked to a quantity x by a relation of the type y = u*exp(v*x)+w.  I know the experimental results (1.186, 0.955, 0.786, 0.648, 0.568, 0.500), to which I cannot assign a precision, associated with the values 1, 2, 3, 4, 5, 6.  I must determine u, v and w.

x=c(1,2,3,4,5,6); y=c(1.186,0.955,0.786,0.648,0.568,0.500)
Plane(0,7, 0,1.3); POINT(x,y, "brown")
yy3 = function(U,V,W,x) U*exp(V*x)+W; Fun = function(x) U0*exp(V0*x)+W0
dif = 1e300
m1=0; M1=10; m2=-100; M2=0; m3=0; M3=1/2; F_RUN3(x,y,1e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
1.091296 -0.3467326 0.3593317  d: 0.0881476985674848 
graph1(Fun,0,7, "red")
m1=0.5; M1=1.5; m2=-1; M2=0; m3=0.2; M3=0.4; F_RUN3(x,y,1e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
1.19855 -0.3138632 0.3137454  d: 0.0113071865924665
...
m1=1.19; M1=1.2; m2=-0.35; M2=-0.3; m3=0.3; M3=0.35; F_RUN3(x,y,1e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
1.190008 -0.3145388 0.318526  d: 0.0104177052956292 
graph1(Fun,0,7, "brown")
# y = 1.190*exp(-0.3145*x) + 0.3185

Three
A weight oscillates on a vertical spring. I know that its height h (in cm) when time t (in tenths of a second) varies is h = U·cos(ωt)+V·sin(ωt).  I can measure ω and find that it is 1.03 s−1. Through a filmed shot I find that to t equal, in order, to 1, 3, 5, 7, 9 correspond approximately to the values of h 3, −16, 6, 9, −8. Estimate how much U and V measure (in cm).

t = c(1,   3, 5, 7,  9)
h = c(3, -16, 6, 9, -8)
range(t);range(h)
#  1 9   -16   9
BF=5; HF=3
Plane(0,10, -20,15); POINT(t,h, "brown")
omega = 1.03; yy2 = function(U,V,x) U*cos(omega*x)+V*sin(omega*x)
Fun = function(x) U0*cos(omega*x)+V0*sin(omega*x)
dif=1e300
m1=-100; M1=100; m2=-100; M2=100; F_RUN2(t,h,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 11.89429 -0.4119553  d: 6.64664848751153
m1=-50; M1=50; m2=-10; M2=10; F_RUN2(t,h,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 12.15634 -0.7992009  d: 6.61732609665133
m1=10; M1=15; m2=-1; M2=-0.5; F_RUN2(t,h,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 12.10602 -0.7896324  d: 6.61680941853707
m1=12; M1=13; m2=-0.8; M2=-0.75; F_RUN2(t,h,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 12.10647 -0.7943228  d: 6.61680536155479
m1=12.1; M1=12.11; m2=-0.8; M2=-0.79; F_RUN2(t,h,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 12.10616 -0.7946137  d: 6.61680532529941
m1=12.106; M1=12.107; m2=-0.795; M2=-0.794; F_RUN2(t,h,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 12.10615 -0.794591  d: 6.61680532519233


# h = 12.11*cos(omega*t)-0.795*sin(omega*t)

Four
What is the function x → A·x^3+B·x^2+C·x whose graph best approximates the points (−5.4), (−2, −1), (−1, −3), (1,5), (3,0), (4, −6)?
If we did not require the graph to pass by (0,0) we could find it with regression3. 

x = c(-5,-2,-1,1,3, 4); y = c(4, -1,-3,5,0,-6)
range(x) ;range(y)
# -5  4   -6  5
Plane(-6,5, -7,6); POINT(x,y, "brown")
yy3 = function(U,V,W,x) U*x^3+V*x^2+W*x; Fun = function(x) U0*x^3+V0*x^2+W0*x 
# In order not to proceed at random I try to understand how (more or less)
# the polynomial must be done (from the points I understand that it has a
# negative directive coefficient and a parameter of degree 3 not too big)
dif = 1e300
m1 = -10; M1 = 0; m2 = -1e3; M2 = 1e3; m3 = -1e3; M3 = 1e3; F_RUN3(x,y,2e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# -1.346109 -3.806864 12.93589  d: 100.149417394534 
graph1(Fun, -7,6, "red")
m1 = -10; M1 = 0; m2 = -1e3; M2 = 1e3; m3 = -1e3; M3 = 1e3; F_RUN3(x,y,2e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# -0.3605754 -1.075944 1.372042  d: 33.967719669027
graph1(Fun, -7,6, "seagreen")
m1 = -2; M1 = 0; m2 = -5; M2 = 0; m3 = 0; M3 = 10; F_RUN3(x,y,2e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# -0.1448519 -0.207785 1.950992  d: 4.62038929996649
# ...
m1 = -0.155; M1 = -0.152; m2 = -0.22; M2 = -0.19; m3 = 1.8; M3 = 2.1; F_RUN3(x,y,2e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# -0.1539919 -0.2125809 1.928769  d: 4.39882991977069 
graph1(Fun,-7,6, "brown")
# y = -0.1540*x^3 - 0.2126*x^2 + 1.929*x

Five
Find the function of the type x → A+B·x+C·sin(π·x/2) whose graph best approximates the points (0,−1), (1,1), (2,1), (3,2), (4,3).

x = c(0,1,2,3,4); y = c(-1,1,1,2,3)
Plane(-1,5, -2,4); POINT(x,y, "brown")
yy3 = function(U,V,W,x) U+V*x+W*sin(pi/2*x); Fun = function(x) U0+V0*x+W0*sin(pi/2*x)
dif = 1e300
m1=-5; M1 = 0; m2 = -1e2; M2 = 1e2; m3 = -1e2; M3 = 1e2; F_RUN3(x,y,2e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# -0.6856656 0.9383885 0.8853998  d: 0.854635046012712
graph1(Fun, -2,6, "red")
m1=-2; M1 = 0; m2 = -10; M2 = 10; m3 = -10; M3 = 10; F_RUN3(x,y,2e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# -0.8583343 0.9894515 0.501509  d: 0.576825403584979 
m1=-1.5; M1 = 0; m2 = 0; M2 = 2; m3 = 0; M3 = 1; F_RUN3(x,y,2e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# -0.8161146 1.00398 0.4829206  d: 0.548950093971122
# ...
m1=-0.82; M1 = -0.78; m2 = 0.95; M2 = 1.05; m3 = 0.45; M3 = 0.55; F_RUN3(x,y,2e5); cat(U0,V0,W0," d:",NUMBER(dif^0.5),"\n")
# -0.7999081 1.00002 0.5004336  d: 0.547722952719513
graph1(Fun, -2,6, "brown")
# y = -0.7999+x+0.50004*sin(pi/2*x)  o  y = -0.80+x+sin(pi/2*x)/2

Six
Find the function of the type x → A·xB whose graph best approximates the points (0.9,1.6), (1.4,2.4), (2,3.3), (3.2,4.4), (4,5).

x = c(0.9,1.4,2,3.2,4); y = c(1.6,2.4,3.3,4.4,5)
Plane(0,5, 0,6)
POINT(x,y, "brown")
yy2 = function(U,V,x) U*x^V; Fun = function(x) U0*x^V0
dif = 1e300
m1 = 0; M1 = 10; m2 = 0; M2 = 10; F_RUN2(x,y,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 1.929549 0.7078059  d: 0.28739501222601 
m1 = 0; M1 = 5; m2 = 0; M2 = 5; F_RUN2(x,y,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 1.909829 0.7070296  d: 0.272840336883223 
m1 = 1.5; M1 = 2.5; m2 = 0.5; M2 = 1; F_RUN2(x,y,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 1.900816 0.711557  d: 0.272237173340978 
m1 = 1.7; M1 = 2.1; m2 = 0.6; M2 = 0.8; F_RUN2(x,y,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 1.898605 0.7124717  d: 0.27221748879405 
m1 = 1.85; M1 = 1.95; m2 = 0.66; M2 = 0.75; F_RUN2(x,y,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 1.898828 0.7123531  d: 0.272217042392632 
m1 = 1.88; M1 = 1.91; m2 = 0.69; M2 = 0.73; F_RUN2(x,y,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# 1.898733 0.7123597  d: 0.272216960655014 
graph1(Fun, 0,5, "brown")
# y = 1.899*x^0.7124  o  y = 1.9*x^0.712
   

Seven
To study the development of myelin as age increases, the total mass of myelin in the brains of mice is measured. The following results are obtained (a: age in days, m: myelin mass in milligrams). The data is taken from an article by D. Agrawal et al. appeared in the "Journal of Structural Biology" in 2009.

a = c(12,16,  17, 22, 57, 135,180,410)
m = c( 0, 3,12.5, 23, 42,  56, 63, 90)

We find the equation of a curve that best fits the data and we use it to determine the mass of myelin at the age of 300 days. I sense that data can be approximated with a logarithmic function.

 
a = c(12,16,  17, 22, 57, 135,180,410)
m = c( 0, 3,12.5, 23, 42,  56, 63, 90)
BF=4; HF=3
Plane(0,500, 0,100); POINT(a,m, "brown")
abovex("days"); abovey("mg")
yy2 = function(U,V,x) U+V*log(x); Fun = function(x) U0+V0*log(x)
# From the data I understand that U must be negative and V must be positive.
dif = 1e300
m1=-1000;M1=0; m2=0;M2=1000; F_RUN2(a,m,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# -57.35715 24.06397  d: 11.7913499114138 
m1=-60;M1=55; m2=20;M2=30; F_RUN2(a,m,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# -57.56779 23.93018  d: 11.6127082100333 
m1=-58;M1=-57; m2=23;M2=25; F_RUN2(a,m,1e5); cat(U0,V0," d:",NUMBER(dif^0.5),"\n")
# -57.81912 23.99737  d: 11.610175340702
# I can stop here. Wanting, I could go on, coming to:
# -57.82579 23.99891  d: 11.6101739459808
# I take:
M = function(t) 24*log(t)-58; graph1(M, 11,500, "seagreen")


M(300)
# 78.89078
POINT(300,M(300), "red")

I can assume that the myelin mass after 300 days is about 79 mg.
 

Seven
We can also use the program to find a curve of a certain type that passes exactly through some known points. An example.

If I know that a function:
F = function(x) exp(-a1*x)*(a2*cos(a3*x)+a4*cos(a5*x))
passes through the points (q, p):
q = c(  0.5,   1.5,  2.5,    3,    4,   5.3,    7,     8,   9.6)
p = c(1.572,-2.026,0.277,0.808,0.193,-0.251,-0.047,-0.193,0.208)
I can estimate the intervals in which the parameters a1, …, a5 fall in the following way (see m1,M1, …, m5,M5).
I choose a high number of steps (2e6). The program will take some time to do all the tests. 

q=c(  0.5,   1.5,  2.5,    3,    4,   5.3,    7,     8,   9.6)
p=c(1.572,-2.026,0.277,0.808,0.193,-0.251,-0.047,-0.193,0.208)
BF=4; HF=3
Plane(-1,10,-3,5)
POINT(q,p, "blue")
yy5 = function(U,V,W,Z,J,x) exp(-U*x)*(V*cos(W*x)+Z*cos(J*x))
Fun = function(x) exp(-U0*x)*(V0*cos(W0*x)+Z0*cos(J0*x))
dif=1e300
m1=0.1;M1=0.5; m2=0;M2=10; m3=0;M3=10; m4=0;M4=10; m5=0;M5=10; F_RUN5(q,p,2e6); cat(U0,V0,W0,Z0,J0," d:",NUMBER(dif^0.5),"\n")
# 0.3454759 0.7815283 9.458656 3.519863 1.945944  d: 0.223866511089548
# Especially at the beginning, I change the extremes of one parameter at a time.
m1=0.2;M1=0.5; m2=0;M2=10; m3=0;M3=12; m4=-0;M4=10; m5=0;M5=10; F_RUN5(q,p,2e6); cat(U0,V0,W0,Z0,J0," d:",NUMBER(dif^0.5),"\n")
# 0.3053948 3.578398 2.012608 1.453747 3.2221  d: 0.182332534262799
m1=0.2;M1=0.4; m2=0.5;M2=6; m3=1;M3=10; m4=-0;M4=4; m5=1;M5=5; F_RUN5(q,p,2e6); cat(U0,V0,W0,Z0,J0," d:",NUMBER(dif^0.5),"\n")
# 0.3275804 3.650718 1.995526 1.453098 3.315538  d: 0.155574765273427 
m1=0.3;M1=0.4; m2=2;M2=4; m3=1.5;M3=2.5; m4=0.5;M4=3; m5=2;M5=5; F_RUN5(q,p,2e6); cat(U0,V0,W0,Z0,J0," d:",NUMBER(dif^0.5),"\n")
# 0.3289414 3.475897 1.993806 1.118655 3.219718  d: 0.0738620087477505
m1=0.3;M1=0.35; m2=3;M2=4; m3=1.9;M3=2.01; m4=1;M4=2; m5=2;M5=4; F_RUN5(q,p,2e6); cat(U0,V0,W0,Z0,J0," d:",NUMBER(dif^0.5),"\n")
# 0.3272992 3.602388 1.999739 1.410829 3.25894  d: 0.0180959724648653
# I can stop here
F = function(x) exp(-a1*x)*(a2*cos(a3*x)+a4*cos(a5*x))
a1=0.33; a2=3.60; a3=2; a4=1.41; a5=3.26
graph2(F,-1,10, "seagreen")