Sunspots formula I used the package 'EUREQA', version Formulize, to analyse the monthly smoothed sunspot timeseries from 1750 till 2010.
It gives me a simple formula, with 8 coefficients, that match data with a correlation coefficient of 0.99214404.
In order to obtain a formula only dependent on past values, to ease the projection into some future, I used dS as a one month delayed difference and asked for a formula of the form: S = f(dS, dS^2 , dS^3).
The best formula obtained is:
   SunSpots = 27.36486757 + dS
    + 27.36486757*(sma(dS, 110))
    + 21.73018064*(sma(dS, 73))
    + 10.03456042*(sma(dS, 34))
    + 2.693356275*(sma((dS^2), 131))
    + 2.189454695*(sma((dS^2), 403))
    + (27.36486757*(sma(dS, 110))*(sma(dS, 110)) - 21.73018064)/(sma((dS^2), 131))

sma(var,length) is the simple moving average.
It reveals a cyclic dependence on 34,73,110,131 and 403 months and a cross between 110 and 131 months (responsable for a long cycle).
The correlation coefficient is quite high, and I expected that the formula could give me some values into the short future. In the last months of the timeseries the error became higher and the formula is quickly divergent in the future.


*

*What can cause such a discrepancy?

*What kind of physical model could give those components in the formula?


Any ideas ?

 A: You've got about 24 sunspot cycles in your raw data; and you've fitted a model with at least 16 coefficients. (there are eight figures that are given to ca. 10 s.f.; and there are eight figures for number of months being averaged. One might make the case of there being additional coefficients too: 1.0 for dS, and the exponents)
So of course your out-of-sample data is going to be well outside the model: all you've done is curve-fitted without any causal basis, and over-specified the model. Fast, large divergence is  exactly what one would expect in the circumstances.
NASA has some helpful thoughts on how to predict sunspot cycles, including references to papers. They report that predicting in-cycle behaviour is "fairly reliable" once the cycle is at least three years in; they use the length and size of one cycle, and the amount of activity at the minimum, to estimate the next cycle. To this, they add data from measurements of changes in Earth's magnetic field, because these changes are caused by solar storms.
I'd start with NASA's modelling, because they have some skin in the game: In their words (from the above link:

"Planning for satellite orbits and space missions often require
  knowledge of solar activity levels years in advance."

A: It is proverbial that overfitting, by using an unparsimonious model, leads to poorer out-of-sample predictions than a more parsimonious model with a worse fit.  Chatfield puts this less strongly: « Although statistics, like the AIC (Akaike's Information Criterion) and BIC (Bayesian Information Criterion), penalize more complex models, the reader should realize that there is still a danger that fitting many models to the same data may give a spuriously complex model that appears to give a good fit, but which nevertheless gives poor out-of-sample predictions.»  The Analysis of Time Series, 6th ed., p. 265.
In the financial world, the mutual funds with the highest returns one year almost always perform more poorly than average the next year.
No one knows why, but this is why you should hire a professional instead of using a software package yourself.
A: My short answer to 

1 - What can cause such a discrepancy ?

Deterministic chaos

2 - What kind of physical model could give that components in the formula?

Deterministic, in  chaos theory, means that there exist dynamical equations that govern the system under study, which individually are absolutely predictive. Acting together on the system they induce a behavior that  to start with seems chaotic in the everyday sense of random, but after studying the results one can find a generalized predictive behavior that cannot be relied on for predicting the outcome of the next time step in an analytic solution. 
