I am fitting a linear regression in R now, and I find that there is one parameter showing linear relationship before certain point and cubic relationship after that point.
What I did is to separate the whole dataset with reference to that point and fit the model with parameter value and cubic value respectively. However, that point is determined manually.
Does anyone here know how to determine that point with formal algorithm?
While you don't specify, I assume you mean the model to be both continuous and smooth at the join.
Such a model might be called a semi-natural cubic spline (natural on one side, ordinary on the other). You can build one by having a linear predictor and a cubic-spline basis function of the form $(x-k)_+^3$ where $k$ is the knot location, and $(.)_+$ means that when the argument is negative the output is 0, while when it's non-negative, the argument is returned unchanged.
e.g. see section 5.2 of Hastie, Tibshirani and Friedman's "Elements of Statistical Learning" 2nd ed (10th printing is downloadable from at least two of the author's web pages) -- see the link here (scroll down if you don't see it).
Your problem reduces to one of automated knot selection (sometimes called knot placement). Hastie et al do discuss this knot selection issue (and problems with it, along with some alternatives to it), that's probably a good place to start. There are many papers relating to this issue (one example is here, though it may not be a suitable reference for you at present).
More prosaically, one such knot can be placed by trying a variety of points; it's equivalent to selecting between the possible candidate locations based on some criterion (such as AIC). However, this has the same problem that model selection more generally has (for which, again, Hastie et al is probably a good place to start, or you could read say chapter 4 of Frank Harrell's Regression Modeling Strategies).
You might want to investigate approaches like cross-validation.
Of course you could always include the sought point into your optimization problem, which then, however, will become non-linear in general and thus more difficult to solve than the normal linear regression problem.
I think one standard approach is to assign a set of points before the calculation and then use piecewiese polynomials in the arising regions (or splines if you need continuity at the boundaries) . Often, the nodes are determined by all or a subset of the data points.
Alternatively, also tree schmemes like CART (classification and regression trees) could be used (with a generalized linear model in the leafs). The approach the optimization problem mentioned before by a greedy algorithm.
