1. A color matching method for a digital image comprising: a step of extracting signals of a plurality of main colors in the digital image as input signals; a step of transformation for generating a 3×3 color matrix from the signals of the plurality of main colors and transforming the signals of the plurality of main colors into a 3×1 input matrix formed by the intensity of illumination, the color hue, and the color saturation of each of the color signals, which are referred to as Lab coordinates, and the color matrix and the input matrix are given as follows: the ⁢ ⁢ color ⁢ ⁢ matrix ⁢ : ⁢ [ a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 ] the ⁢ ⁢ input ⁢ ⁢ matrix ⁢ : ⁢ [ L i ⁢ ⁢ n a i ⁢ ⁢ n b i ⁢ ⁢ n ] wherein each of the elements a 11 to a 33 in said color matrix is a real number, and L in , a in , and b in are respectively representative of the transformed Lab coordinates of each of the color signals; and a step of color-difference calculation for recursively calculating an output value and color difference of the plurality of the color signals by performing an iteration operation based on Equation (1) composed of said color matrix, a weighting matrix and a fine-tuning matrix, until said color difference is smaller than a predetermined threshold value, said color difference is defined as: Δ ⁢ ⁢ E ⁢ ⁢ ( color ⁢ ⁢ difference ) = ( L out - L Ref ) 2 + ( a out - a Ref ) + ( b out - b Ref ) 2 wherein L out , a out , and b out are representative of the obtained Lab coordinates for one of said colors, L Ref , a Ref , and b Ref are representative of reference values of the obtained Lab coordinates for said one color and the Equation (1) is given as follows: [ L out a out b out ] = ( [ a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 ] + [ w 11 w 12 w 13 w 21 w 22 w 23 w 31 w 32 w 33 ] ⁡ [ k 11 k 12 k 13 k 21 k 22 k 23 k 31 k 32 k 33 ] ) × [ L in a in b in ] Equation ⁢ ⁢ ( 1 ) wherein [ L out a out b out ] is the outputted Lab-coordinate matrix of said one color, [ w 11 w 12 w 13 w 21 w 22 w 23 w 31 w 32 w 33 ] is said weighting matrix, and [ k 11 k 12 k 13 k 21 k 22 k 23 k 31 k 32 k 33 ] is said fine-tuning matrix wherein each of the elements in said fine-tuning matrix is defined as follows: 
 k ij =P ij +D ij +I ij 
 wherein 
 P ij =Ep=K p ×ΔE ( t ) 
 D ij =Ed=K d ×(Δ E ( t )−Δ E ( t− 1)) 
 I ij =Ei=K i ∫ −∞ ΔE ( t ) dt wherein K p , K d , and K i are representative of the coefficient of proportional control, the coefficient of derivative control, and the coefficient of integral control.